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ABSTRACT 

The  unsteady,  sonic,  small  disturbance  potential  equation 
was  approximated  by  reduction  to  parabolic  form  and  solved  by 
Laplace  and  Fourier  Transform  and  Collocation  techniques, 
including  thickness  effects  through  the  Oswatitsch  parameter. 
The  Transform  techniques  produced  reliable  solutions  and  the 
Collocation  solution  was  usable  with  qualifications,  for  some 
parameter  values,  stemming  from  the  choice  of  approximating 
series.   Stability  boundaries  for  an  unstaggered  cascade 
oscillating  in  sonic  flow  were  calculated  with  the  Fourier 
Transform  technique  revealing  that  blade  interference  is 
strongly  destabilizing.   The  Collocation  technique  was  suc- 
cessfully applied  to  the  two-dimensional,  unsteady,  porous 
walled,  sonic  wind  tunnel  interference  problem,  heretofore 
unapproached  with  the  inclusion  of  thickness  effects. 
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I.   INTRODUCTION 

Nonsteady,  interfering,  transonic  flows  are  of  increasing 
importance  to  aeroelasticians.   Recent  design  concepts  for 
transonic  transports,  high  maneuverability  fighters,  biplanar 
space  shuttle  configurations,  and  so  called  transonic  axial 
flow  compressors,  to  name  a  few,  involve  such  flows  when  aero- 
elastic  problems  are  considered.   Ironically,  while  sub/super- 
sonic flutter  problems  have  been  extensively  studied,  the 
transonic  problem  has  been  circumvented  by  either  rapid 
acceleration  through  the  regime  or  operation  avoiding  it.   It 
is  now  apparent  that  the  question  of  the  effect  of  transonic 
flow,  including  possible  adjacent  flow  field  interference, 
on  flutter  must  be  confronted  in  a  practical  way. 

Of  course,  even  without  interference,  nonsteady  transonic 
flows  challenge  the  aeroelastician  concerned  with  aerodynamic 
damping  of  structural  vibrations.   In  that  regime,  pressure 
perturbations  generated  by  a  body  and  propagated  at  the  local 
speed  of  sound  tend  to  accumulate  near  the  body  itself.   Thus, 
the  magnitude  of  aerodynamic  forces  is  increased  and,  perhaps 
more  importantly,  a  large  phase  lag  can  exist  between  the 
forces  and  the  motion  producing  them.   It  is  this  phase  lag 
that  can  lead  to  negative  aerodynamic  damping  and,  potentially, 
dangerous  flutter. 

It  is  not  enough  to  analytically  study  the  flutter  charac- 
ter of  particular  configurations.   Wind  tunnel  experiments 
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must  be  conducted  to  verify  analytical  results  before  a 
configuration  can  be  committed  to  flight  test.   However, 
flow  perturbations  from  a  wind  tunnel  test  model  may  reflect 
from  the  tunnel  walls  back  onto  the  model  requiring  a  correc- 
tion to  test  results  due  to  interference.   In  steady  transonic 
flow,  such  wall  interference  has  been  thoroughly  examined  and 
porous  or  slotted  wall  tunnels  have  been  utilized  for  some 
time  to  reduce  its  effect.   In  nonsteady  flow,  the  porous 
wall  boundary  condition  is  a  difficult  one  for  classical 
analytical  approaches  to  handle.   Only  Drake  [Ref.  6]  has  in- 
cluded wall  porosity  in  a  theoretical  analysis  of  a  thin  air- 
foil oscillating  in  a  transonic  wind  tunnel  and  his  results 
are  by  no  means  complete.   Thus,  there  is  a  great  need  for 
solutions  that  can  yield  the  magnitude  of  unsteady  wall 
interference  and  the  effect  of  wall  porosity  upon  it. 

Blade  flutter  in  transonic  axial  compressors  is  a  complex 
problem  involving,  among  other  things,  interference  of  highly 
perturbed  transonic  flow  fields.   Historically,  sub  and  super- 
sonic axial  compressor  flows  have  been  simulated  by  a  series 
of  airfoils  (or  blades)  placed  in  close  proximity  (a  cascade) 
so  that  the  two  dimensional  flow  models  that  in  a  rotating 
compressor  stage.   For  the  transonic  compressor  stage  the 
mixed  flow  character  of  the  real  flow  makes  the  two-dimen- 
sional approximation  questionable,  e.g.,  McCune  [Ref.  16]. 
It  must  be  recognized  that  accurate  modeling  of  such  a  complex 
flow  field  is  well  beyond  the  state  of  the  art.   However,  the 
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cascade  approximation  may  yield  important  information  relative 
to  the  trends  of  blade  flutter. 

Nonsteady  transonic  flow  analysis  begins  with  the  deriva- 
tion of  the  nonlinear  small  perturbation  potential  equation 
which  Lin,  Reissner,  and  Tsien  [Ref.  15]  and  Landahl  [Ref. 
13]  have  shown  can  be  linearized  for  sufficiently  high  fre- 
quency.  Since  aeroelastic  analyses  are,  generally,  concerned 
with  the  stability  of  small  amplitude  oscillations,  the  solu- 
tions to  the  nonlinear  potential  equations  shall  be  written 
as  the  sum  of  a  steady  component  and  a  much  smaller  nonsteady 
component.   The  result  is  two  equations,  the  usual  nonlinear 
steady  flow  equation  and  a  linear  nonsteady  flow  equation 
containing  variable  coefficients  depending  upon  the  steady 
flow. 

Transonic  cascade  flow  for  thin  blades  and  in-phase  oscil- 
lations was  approached  by  Hamamoto  [Ref.  10]  using  Fourier 
Transform  techniques  on  the  linearized  potential  equation. 
Gorelov  [Ref.  8]  solved  the  linear  nonsteady  equation  for 
mildly  supersonic  flow  past  thin  blades  with  a  superposition 
of  isolated  blade  and  interference  potentials  represented  as 
trigonometric  series.   See  also  Samoylovich  [Ref.  22]. 
Elder  [Ref.  7]  revised  Gorelov 's  Collocation  technique  by 
using  a  power  series  and  solved  for  the  steady  sonic  flow 
cascade  potential  and  the  coefficients  of  the  series  for  the 
nonsteady  flow,  out-of-phase  oscillation  case.   Elder  made  use 
of  the  Oswatitsch  [Ref.  20]  parameter  (used  by  Oswatitsch  for 


steady  flows),  as  extended  by  Teipel  [Ref.  30]  and  Hosokawa 
[Ref.  11]  for  nonsteady  flows,  to  allow  for  blades  of  nonzero 
thickness  (nonuniform  steady  potential  flow). 

Sonic  wind  tunnel  interference  in  steady  flow  was  solved 
by  Sandeman  [Ref.  23]  using  Spreiter's  local  linearization 
technique,  this  after  Guderley  [Ref.  9]  and  Spreiter  et  al 
[Ref.  27]  showed  that  the  choked  tunnel  results  can  be  cor- 
rected so  as  to  simulate  unbounded  flows.   Murasaki  [Ref.  18] 
solved  steady  transonic  tunnel  flows  for  porous  or  slotted 
walls.   Miles  [Ref.  17]  and  Drake  [Ref.  5J  solved  the  super- 
sonic nonsteady  solid  and  porous  wall  cases,  respectively. 
Drake  [Ref.  6],  then,  utilized  Laplace  Transform  techniques 
and  the  linear  potential  equation  for  a  thin  airfoil  between 
porous  walls  for  transonic  flows.   Recently,  Savkar  [Ref.  2*4] 
has  Laplaced  transformed  the  linearized  equation  and  studied 
airfoil  stability  in  a  solid  walled  wind  tunnel. 

Wind  tunnel  interference  can  be  approached,  according  to 
Sears  [Ref.  25],  with  suitably  located  sensors  within  the 
test  section  (not  on  the  model)  and  variable  wall  shape  or 
porosity  and/or  plenum  pressure  behind  the  wall.   The  sensors 
would  modify  the  tunnel  parameters  until  unbounded  flow  con- 
ditions were  present  at  them.   Sears  argued  that  unbounded 
flow  conditions  are  thus  obtainable  upon  the  model. 

This  work  began  with  the  object  of  applying  a  Collocation 
technique  to  the  porous  wall  wind  tunnel  problem  with  the  con- 
viction that  that  method  could  suitably  handle  the  classically 
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difficult  porous  wall  boundary  condition.   While  proceeding 
toward  that  goal,  an  efficient  approximate  solution  was  de- 
veloped for  studying  the  aerodynamic  stability  of  the  sonic 
oscillating  cascade  and  stability  boundaries  were  calculated 
for  the  cascade  as  well  as  for  a  biplanar  configuration. 

To  accomplish  the  initial  objective  of  this  work,  the 
nonsteady  potential  equation  was  reduced  with  the  aforemen- 
tioned Oswatitsch  parameter  to  approximate  the  variable 
coefficient  therein.   The  resulting  reduced  equation  was 
solved  by  a  Collocation  technique  for  the  cascade  problem 
following  Elder  [Ref.  71.   For  numerical  checks  upon  the 
Collocation  technique,  the  equation  was  Laplace  transformed 
extending  the  steady  flow  work  of  Sandeman  [Ref.  23]  and 
Fourier  transformed  in  a  manner  analogous  to  Hamamato  [Ref. 
10].   These  extensions  were  made  flexible  enough  to  handle 
blade  thickness  variations  (through  the  Oswatitsch  parameter) 
and  variable  interblade  phase  angle. 

The  Elder  power  series  initially  utilized  in  the  Collo- 
cation solutions  to  the  cascade  problem  was  found  to  have 
numerical  difficulties  that  increased  computation  time  re- 
quirements and  reduced  accuracy.   Thus,  a  'step-wise-linear' 
series  was  introduced  and  the  resulting  Collocation  method 
solution  to  the  cascade  problem  agreed  well  enough  with  the 
other  two  solutions  to  allow  some  confidence  in  its  subsequent 
adaptation  to  the  porous  wall  wind  tunnel  problem. 
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The  Collocation  solution  was  found  to  be  quite  flexible. 
It  can  be  altered  to  allow  streamwise  variation  of  the  Os- 
watitsch  parameter,  in  a  manner  reminiscent  of  Spreiter's 
local  linearization,  to  improve  accuracy.   Further,  it  allows 
a  variation  in  wall  porosity  for  possible  simulation  of  the 
technique  mentioned  by  Sears. 
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II.   UNSTEADY  TRANSONIC  PLOW  THEORY 

A.   EQUATIONS  FOR  GENERAL  OSCILLATORY  MOTION 

The  governing  equation  considered  below  arises  from  the 
assumption  that  the  flows  of  interest  are  of  frictionless 
perfect  gas.   Thence,  the  equation  of  state  and  Eulerian 
equations  of  motion  are  invoked  and  simplified,  assuming  the 
absence  of  heat  conduction  and  strong  shock  processes.   Fur- 
ther, a  velocity  potential  is  introduced  which  reduces  the 
number  of  dependent  variables  at  the  imposition  of  the  re- 
quirement of  flow  irrotationality.   Finally,  a  small  distur- 
bance analysis,  including  the  unique  properties  of  transonic 
flows  (see  Refs.  12  and  13),  yields,  in  two  dimensions 


<J>TT   +    2(J>XT   +    [1    -   -^  +    (Y+l)4>x]cj>xx   -   -2<J>YY    =    0.         (II-l) 


This  equation  is  applicable,  approximately,  to  flows  containing 
weak  shocks,  incidentally  (see  Ref.  13). 

In  the  study  of  unsteady  aerodynamics,  it  is  customary 
and  convenient  to  expand  the  solution  into  steady  and 
oscillatory  nonsteady  components.   Thus, 


(J>(X,Y,T)  =  4>U,Y)  +  R.P.[>(X,Y)eiKT],         (H-2] 


where  K  is  the  reduced  frequency. 
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Substitution  in  equation  (II-l)  yields  a  coupled  partial 
differential  equation, 


(1-M2)tJ>xx  +  4)yY  -  M2  (y+1)  T^x^xx  +  *x^xx  +<^x^x] 

+  M2K2^  -  2iM2Ki^Y  =  0.  (H-3) 


This  equation  is  nonlinear,  of  mixed  elliptic-hyperbolic 
type,  and  contains  variable  coefficients.   However,  because 
flutter  analysis  is,  generally,  concerned  with  the  stability 
of  small  perturbations  about  a  steady  flow,  it  is  appropriate 
to  assume  that  the  oscillatory  nonsteady  component  is  small 
compared  to  the  steady  flow  potential.   Thus,  the  governing 
equation  becomes 


(1-M2HXX  +    ipyY   -   M2(Y+l)[<t>x^xx  +    <f>xx^xl 

-  2iM2K^v  +  M2K2ip  =  0,  (H-4) 


with  the  usual  nonlinear  steady  flow  equation  also  resulting. 
Equation  (II— 4) ,  though  linear,  retains  variable  coefficients 
dependent  upon  the  steady  flow  solution,  the  nature  of  which 
(subsonic  or  supersonic)  determines  its  type  (elliptic  or 
hyperbolic) . 

B.   REDUCTION  OF  THE  GOVERNING  EQUATION 

First,  the  reference  Mach  number  shall  be  unity,  reducing 
equation  (II-lJ)  to 
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^YY  ~  ("Y+DUx^x  +  ^XX^X^  "  2iK%  +  r2^'  =  °"         (H-5) 


Second,  the  approach  taken  by  Teipel  [Ref.  30]  and 
Hosokawa  [Ref.  11]  is  utilized.   That  is,  the  nonconstant 
coefficients  on  i^„„  and  ty„   are  approximated  as  follows. 
First,  the  steady  flow  velocity  perturbation  is  assumed  to 
be  small  in  the  X  direction,  implying  that  the  airfoil  or 
blade  is  of  small  thickness, 


4>x  1   0.  (II-6) 


Second,  the  steady  flow  fluid  acceleration  over  the  airfoil 
or  blade  is  assumed  to  be  a  constant  (parabolic  arc  airfoils 
exhibit  nearly  this  characteristic),  which  allows  the  intro- 
duction of  the  Oswatitsch  parameter 


X  =  (Y+D<j>xx  =  constant  >   °-  (H-7) 


The  Oswatitsch  parameter  was  originated  by  Oswatitsch  and 
Keune  [Ref.  20]  for  steady  flow,  incidentally.  Thus,  the 
governing  equation  has,  finally,  become 


\\>YY   -    (X+2iK)i|;x  +  K24>   =  0.  (H-8) 


Note  that  this  equation  is  of  the  parabolic  type. 
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C.   BOUNDARY  CONDITIONS 
1.   Flow  Tangency 

The  usual  surface  boundary  condition  is 


H  =  0,   on   F(X,Y,T),  (II-9) 


DT 


where  F(X,Y,T)  describes  the  surface  of  a  body  moving  as  a 
function  of  time.   The  boundary  condition  requires  flow  tan- 
gency, that  is  the  total  derivative  of  F  is  zero.   Then,  con- 
sidering a  thin  airfoil  restricted  to  vanishingly  small  ampli- 
tude out  of  plane  oscillations  (recall  the  small  nonsteady 
potential  argument) ,  its  upper  and  lower  surfaces  can  be 
described  by 


F  (X,Y,T)  =  Y  -  H  (X,T) , 
u   5  '  u   '   ' 


F£(X,Y,T)  =  Y  -  H£(X,T), 


(11-10) 


where 


H  (X,T)  =  xh  (X)  +  R.P.Ch  (X)eiKT]  (11-11) 

u  '       u  u 


With  the  definition  of  4>(X,Y,T)  from  equation  (II-2) 
and  H(X,Y,T)  from  above,  the  linearized  flow  tangency 
condition  becomes 
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♦  l 


.(X,0±)  =  ^P-  +  iKh(X).  (H-12) 


2 .   Zero  Upstream  Disturbance 

The  governing  equation  has  become  parabolic  as  shown 

by  equation  (II-8),  requiring,  upstream  of  the  leading  edges, 

*(X,Y)  =  0.  (H-13) 
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III.   FORMULATION  OF  SOLUTIONS 

A.   SOLUTION  THROUGH  THE  LAPLACE  TRANSFORMATION 
1.   Cascade  Geometry  and  Boundary  Conditions 


u. 


k 


-2b- 
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---  hn 


Fig.  1.   Unstaggered  Cascade. 

Consider  a  pair  of  the  parabolic  arc  airfoils  or 
blades  in  an  infinite  unstaggered  cascade  (row)  immersed 
within  a  transonic  flow  field.   The  blade  chord  is  2b,  inter- 
blade  spacing  is  P,  the  lower  blade  is  the  reference  blade 
(variables  associated  with  it  are  subscripted,  zero),  and 
the  upper  is  referred  to  as  the  adjacent  blade  (subscript, 
one).   The  coordinate  reference  is  fixed  to  the  reference 
blade  midchord  for  consistency  with  Teipel  [Ref.  30]  with 
whose  work  comparisons  were  made.   Pitch  and  plunge  ocsilla- 
tions  can  be  described,  pitch  axis  at  a,  magnitude  of  pitch 
aQ,  of  plunge  h  .   The  adjacent  blade  performs  a  harmonic 
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pitch/plunge  oscillation  arbitrarily  phased  with  the  oscilla- 
tion of  the  reference  blade. 

Two  n  coordinates  are  specified,  with  n0  and  nn 
related  such  that 


n1   =  n0  +  P.  (in-i) 


The  following  nondimensional  variables  are  defined: 

h 

X  =  1T  • 


*™  -  tt  -  m  =  0,1, 


(III-2) 


m    b 


P 
b 


A0   b  ' 


hl  "  "F 


Note  that  the  linear  nondiniensionalizing  parameter 
will  be  the  blade  semichord  for  all  variables  in  this  analysis 
Thus,  for  example,  the  reduced  frequency 


K  -  {&  .  (IH-3) 

1 
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From  considerations  in  Section  II,  the  boundary  con- 
ditions can  be  expressed,  for  the  reference  and  adjacent 
blades,  as 


8h 
iPY(X,YQ  =  0)  =  vQ(X)  =  -^  +  iKhu  ,  (111-4) 


and 


i|;y(X,Y0  =  p)  =  Vl(X)  =  eiav0(X)  ,  (III-5) 


where  a  is  the  interblade  phase  angle. 

Then,  for  pitch  and  plunge  of  small  magnitude 


vQ(X)  =  -  a0[l+iK(X-XQ]  +  iKh_.  (III-6) 


For  computations  for  the  steady  flow  case,  a 

p 
parabolic  arc  airfoil  is  assumed  (Y  =  t(1-X  ),  see  Ref.  30) 

and 


vQ(X)  =  -  2tX.  (HI-7) 


2 .   Solution  of  the  Oscillatory  Sonic  Flow  Equation 
The  governing  equation  from  Section  II  is 

*YY  "  Ql|;X  +   K^  =  °>  (IH-8) 
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where 


Q  =  A  +  21K.  (HI-9) 

For  utilization  of  the  Laplace  Transformation  define 
a  new  streamwise  coordinate 


X.  «  X  +  1.  (111-10) 


Denoting  the   Laplace   transformation   of  a   function 
with  respect   to   X,    with   an   overbar,    i.e., 

*f  °°  "pXi 

f(p)    -«af(X,)]   -     /  e        1f(X1)    dX,    ,  (111-11) 

1  0  ii 

transform  equation    (II-8)    to 


\pYY   -    \i2ip  =    0    ,  (111-12) 


where 


y2    =    Qp   -    K2.  (111-13) 


The   boundary    conditions  become 


^Y(P,0)    =   vo(p),  (111-14) 
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*y(P,p)  =  elavQ(p).  (IH-15) 


The  solution  to  equation  (111-12)  subject  to 
(111-14  and  15)  is 

^(pjY)  =  l^EL   elgcosh(yY)  -  cosh[y(Y-p)]  (III-16) 

vv^5       y  sinhlyp) 

The  hyperbolic  functions  are  expanded,  as  follows, 
to  allow  a  term  by  term  inversion 


_1_  =  _§ m    2  -yp      e-2nyp  (111-17) 

sinh(yp)    (eyp  _  e~yp)         n=Q 


then 


cosh^Y)  =   "  {  e-y[-y+(2n+l)p]  +  e-y[Y+(2n+l)p]}     (III.l8) 


sinh(yp)    n=Q 


and 


cosh[ 
sinh 


y(Y-pli=      ^    {e-y[2np+2p-Y]    +    e-M[2np+Y]  (IH-19) 

h(yp;  „_nL 


n=0 


Thus, 


^(p>Y)    ,   I     ~    {-    (p)eic[e-y[-Y+2nP+p)    +    e-y (Y+2np+p) 3 

y  n=0      u 

(111-20) 

-v0(p)[e"y(-Y+2nP+2P)    +   e-»«^y+2nP)]}. 
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For  inversion,  recall  the  formula 

r  e-pva/i4xi        ,-w 

n -§- dX,  =  -Jv- ,   R.P.[p]  >  0.      (111-21) 

Utilizing  the  above  formula  and  the  convolution 
theorem,  recognizing  that 

vQ(X1)  =0,   Xx  <  0  ,  (111-22) 

and  reintroducing  X   after   inversion,   yields 


K2 
x     ^j-(X-S)      M  _Q(2HP+P-Yil 

*(X,Y)    =  -i-        /     S z    {vn(S)ela[e   ^      (X"S)        + 

IprQ1  -1  -VX^       n=C      U 

Q(2np+p+Y)2  Q(2np+2p+Y)2  Q(2np^Y)2 

e-T»       rx^sT-3  .  Vo(s)[e-^       OPsT—  +  eT— TX=S7]}dS- 

(111-23) 


Then,    for   the   potential   on  the   reference   blade 
upper   surface, 


2  2 

X      £_(x-S)    oo  _Q(2np+p) 

Y(X,0+)    =  -5-        /      e^ Z    (vr,(S)ela[2e    ^  TtST  3 

-vW  -l      Vxls^     n=0 

Q(2np+2p)2  Q(2np)2 

-   vQ(S)[e  ^      U"S)      "  +   e"4    (X-S}]}d3.  (111-2*0 
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This  equation  was  programmed  for  digital  calculation  (see 
Appendix  A  and  Section  IV-E). 

3.   Modification  for  Sonic  Free  Jet  Problem 

The  sonic  free  jet  is  of  interest  as  an  extreme  case 
of  wind  tunnel  wall  porosity.   That  is,  the  porosity  is  so 
great  that  the  wall  has  vanished  and  the  airfoil  is  in  a 
sonic  jet  constrained  at  its  boundaries  not  by  a  zero 
velocity  condition  but  zero  pressure  differential. 

The  boundary  condition  on  the  reference  blade  remains 
the  same  as  before,  equation  (III-4).   The  adjacent  blade  is 
replaced  by  the  edge  of  the  sonic  jet  upon  which  the  following 
pressure  perturbation  condition  is  imposed 

iK^(X,Y=p)  +  ^x(X,Y  =  p)  =  0.  (IH-25) 

Transforming  and  substituting  as  before,  the  second 
boundary  condition  becomes 

(iK  +  p)^(p,p)  =  0.  (111-26) 

Then,  using  both  conditions,  the  solution  to  the 
transformed  equation  is 

v  (p)   -2yp  uY    p-yY 
f(p,Y)  =  — [£ o  ] 

m   Vo  p)[slnhCy(Y-o)]:i  (111-27) 

u   L   cosh(yp) 
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With 


2e~yL   =  2e"yp   Z  (-l)ne"2nyp  ,        (111-28) 


cosh(pP;    (l+e"2yp)         n=0 


?(P,Y)  =  IS^l  "  (-1)n{e-y(,.Y+2np+p)  _  e-u[Y+2np]}<  (m_29) 
y   n=0 


The  subsequent  inversion  is  similar  to  that  in  the 
previous  section. 

The  result  was  programmed  for  digital  computation  and 
used  as  a  check  case  against  which  the  Collocation  technique 
porous  wall  wind  tunnel  solution  was  compared.   The  Fourier 
Transformation  method  (see  next  section)  was  also  applied  to 
the  free  jet  problem  and,  by  following  a  derivation  similar 
to  that  above,  it  was  found  that  an  interblade  phase  angle 
of  zero  in  the  cascade  problem  correctly  modeled  the  sonic 
free  jet  problem. 

B.   SOLUTION  THROUGH  THE  FOURIER  TRANSFORMATION 
1.   Cascade  Geometry  and  Boundary  Conditions 
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Fig.  2.   Unstaggered  Cascade. 

The  arrangement  of  the  cascade  is  the  same  as  the 
previous  section,  however,  the  coordinate  system  is  changed 
to  simplify  the  derivation.   The  origin  of  the  coordinate 
system  is  now  fixed  to  the  reference  blade  leading  edge. 
Blade  chord  is  c. 

The  following  nondimensional  variables  are  defined 
for  this  phase  of  the  study: 


X  = 


Y  =  a 

1    c  » 


(111-30) 


h  = 


c  3 


1    c  y 


3^ 


V,   -  ^ 

h2  "  T 


All  variables  will  use  c  (vice  b)  for  nondimension- 
alization  as  appropriate,  i.e., 


K-,  -  £p  .  (Hl-31) 

±  1 


The  boundary  conditions  are  the  same  as  those  expressed 
in  the  previous  section,  though  their  representation  is 
somewhat  different. 


if>y(X,Y=0)  =  w(X)  , 


^y(X,Y=h)  =  e10w(X)  ,  (IH-32) 


and 


w(X)  =  [-1-1K  (X-A1)]oQ  +  iK1h2  .       (HI-33) 


2.   Transformation  of  the  Sonic  Flow  Equation 

As  before,  begin  with  the  equation  for  oscillatory 
sonic  flow  from  Section  II  and  transform  it  with  the 
following  Fourier  Transform, 


+  0O 

f(v,Y)  =   /   f(X,Y)eiVX  dX  ,  (111-3*0 
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Obtain 


and 


or 


4>YY  +  P2ip  =  0  ,  (111-35) 


y2  =  -  iAv  +  2K  v  +  K±2 .  (III-36) 


This  equation  is  satisfied  by 


J(v,Y)  =  EieiyY  +  F  e""lyY  .  (IH-37) 


Using  equation  (111-32) 


i/jY(v,h)  =  e   i|>y(y,0)  , 


El[eiyh  -  eia]  =  Fl[e-iyh  -  eia]  , 


.  -la  iyh 

e1  =  f.c1  ~  e  .  e  .  .  ]  .  (in-38) 

1    It  -ia  -iyhJ 

1  -  e   e   H 


Then,  with 


ipy(v,0)  =  w(v) , 
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E,    and   F,    become 


_     .  w(v)      rl   -   e-1^-1"11, 
1   "  T^P  Vyh  -   e"luh        ' 


and 


— /    s      n  -la   iuh 

_  w(v)    rl   -   e        e    M    n 

i "  3 — i^[tt ^r] 

iye  iyh  -iyh 


Finally, 


— /    %       /,      -la  -iyhN    iuY    ,     ,_      -ia   iuh\    -iuY 
^(VjY)    =  w(v?    [Hzg e  e    *     +    (1-e e    Me      V      j 

V        iyela  e^h  -  e"^h 

(111-39) 


Because  the  pressures  on  the  reference  blade  are  of 
interest,  set  Y  =  0.   With 


ia,  -ia                 ia  _ia 
e   +e             .      e   _e 
cos  a  =  p ,       sin  a  - 


and 


2 

^„  o   -  1  -  tan  a  n      _        2  tan  a 

cos  2a  =  p— ,  sin  2a  =  p—  » 

1  +  tan  a  1  +  tan  a 


then 


a  .,      9  yh 

*(v,0)  =  -  w^d-e-1")  t   2(l+e-lg)tan2  "Tj  , 
yela         Han  H? 
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or 


?(v,0)  =  -  Hi^i  [{si!=l)00te  +  (S^±-i)tan^].      (III-itO) 


It  can  be  shown  that  equations  (III-16)  and  (111-40) 
are  equivalent  and  in  this  sense  the  two  solutions  are  not 
independent.   However,  the  inversions  are  substantially 
different.   The  following  stems  from  the  work  of  Hamamoto 
[Ref.  10]  and,  apparently,  arises  from  observations  of  abstract 
complex  variable  theory.   It  differs  from  standard  Fourier 
Transform  inversion  techniques  and  leads  to  a  very  compu- 
tationally convenient  result.   A  further  method  of  inversion, 
utilizing  the  Theorem  of  Residues,  was  considered  and  is 
discussed  in  Section  V-C. 
3.   Inversion 

In  order  to  simplify  the  inversion  of  equation  (111-40), 
SN  and  TN  will  be  introduced  and  evaluated  in  a  subsequent 


section. 


Define   SQ   and  Tfl   as    follows, 


+00 
SQ(X)    =  X       f       han  If  eivX  dX,  (111-41) 

_00 


and 


T    (X)    =     *         /        Icot   ^  eivX  dX.  (111-42) 


_0O 
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Also 


4> 
SN+1U)  =  /  sN(n)  dn,  (Hi-4.3) 


and 


4> 
TN+1U)  =  /  TN(n)  dn.  (Ill-W 


Using  the  convolution  theorem  and  recognizing  that 


w(X)  =  0,   ip(X,Y)  =  0,   X  <  0,  (111-45) 


yields 


1|>(X,0)  =  -  /  w(X-4»)[(S-g=l)T0(*)  +  (e_-+l)S0(4,)]  ^ .  (111-46) 


Integrate  by  parts  to  obtain 


X 
Tp(X,0)  =  -{w(0)[D-T_(X)  +  B-S  (X)]  +   /  w(X-(J>)[D-T1(<j>) 

1  1  Q  1 


+  B-S  (<|»)]  d^}  ,  (111-47) 


where 


B  =  e   "  1,    D  =  e  J-   k  (III-U8) 
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4.   Calculation  of  Lift  and  Moment  Coefficients 

This  approximate  solution  of  the  sonic  unsteady 
cascade  problem  (111-47)  was  found  to  be  quite  efficient  in 
its  usage  of  computer  time.   Thus,  it  was  extended,  as  follows, 
to  calculate  pressure,  lift,  and  moment  coefficients  in 
addition  to  velocity  potential. 
Since 

CP(X,0)  =  -2[i|;Y(X,0)  +  iK.^(X,0)]  ,  (111-49) 


and 


CT  =  -2   /   CD(X,0)  dX  ,  (111-50) 

L      0    P 


then 


1 
CL  =  4   /   Ux(X,0)  +  iK1^(X,0)]  dX 


1  1 

=  4[  /   iMX,0)  dX  +  iK.   /   i|;(X,0)  dX].   (III-5D 

o    A  -1  0 


With  the  previously  defined  ty   and  w(X) ,  and  the 
following  definition 


/   4>(X,0)dX  =  F2(l)  =  -{(-1+1K1A1)[D-T2(1)+B-S2(1)] 


-IK1[D.T3(1)+B-S3(1)]}  ,    (111-52) 
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the  lift  is 


CL  =  40(1,0)  +  JXjFgCl)]. 


For  the  moment  coefficient 


CM  =  2   /  (X-A1)  C  (X,0)  dX 


1  1 

=  -4[iK..   /  Xi|;(X30)dX  +   /  X^v(X,0)dX]  +  fl-CT  .     (III-5M 
1  0  0    A  ■  1    L 


Integrating  by  parts  and  continuing  the  notation 
started  above 


CM  =  -H{-iK1[F2(l)-P3(l)]  +  [ip(l,0)-F2(l)]}  +  A1.CL.  (111-55) 


5 .   Solutions  for  S   and  T«j 

Because  SQ(X)  and  TQ(X)  are  difficult  to  evaluate 
directly,  I(X)  and  J(X)  will  be  defined  and  calculated.   It 
will  be  seen  that  I(X)  and  J(X)  are  S.(X)  and  T-^X),  respec- 
tively, from  which  all  SN  and  T   can  be  obtained. 
Define 


,      +«>  ivX      . 

I(X)  =  J±-       f       S tan  &■   dv.  (111-56) 

2tti        uv       2 

-.00  r 


Since  [Kef.  12] 


kl 


fflM    ttX  ^  4X 

Ian      2      —        2j ?j ^-     j 

m=0    (2m+l)"-X 


with 


ttX   _    uh  v   _    uh  , 

-K-  =  "V  ,        X   =  -^   ,  m  =  n-1, 

2  2'  7T       '  ' 


1     ,.orl     TTh     _  1        " 

—  tan  -s-  -    -  t-      Z 


u  2  h -  n  n      _    0      2    s 

n-l    {#)2  _    (2n-lj= 


2>      "    ^2     '      ? 


then 


,  +00     ivX    °°  . 

I(X)    =   -   p~        /        ^-r-[  £       1 ]    dv 

n-l    (U.}2   _    (2r^l)2    tt_ 

h 


Substitute    for   u,    using  equation    (111-35)    in   the 
bracket 


00  00 

I  E    i =  J_    z i 

2      ,       "  ,        .     ~         ,    "    2  K 


5 


n=l   -iXv   +    2Knv  +   Kn  -      _    0      2        "l  n=l,,    IX    w  % 
1 1              ,2n-1^2    -ft                          (l-oiT~Hv-vinJ 

~ 1 -    (~2_)      JT  2K1  n 


=    [(2n"12    *      -   K12]/[2K1-iA]    ,  (IH-57) 

n  h^  i  l 


so 


IT  oo+oo  IvX 

I(X)    =   -  ~-   (ovl   ,  ,)      E  /        -?- r  dv.  (111-58) 

iirh      2K1-iX     n=1  -eo  vU^TT 
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Now   define 


+«     ivX 
J(X)    =  2¥i       /       ^T  cot     2    dv*  (IH-59) 


.CO 


With    [Ref.    12] 

00 

ttX    cot    ttX   =    1   +    2X2      Z      ~2-^ g-    , 

m=l   X      -   v 

obtain 


1  +«   ivX     2  -, 

27T1   -co  vh      y2        n=2    ,^2         (n-l)V 

h 


With 


2    2 
^   =    [iiHzl)_L_  _    K12]/[2K1   -    IX]    ,  (111-60) 


then 


+oo     ivX 


J(X)    ■  T^h    (2K7TTX)      *      7        '        5&=rCT  dV-  (III'61) 

JL  n~JL  -co  j\ 

The    integral   occurring   in   equations    (III-58   and   6.1) 
is,    from  Campbell   and  Foster    [Ref.    4], 

+~      ivX  0    .         iv   X 

'       71^7-7  dv  -  §^  [e     n     -  1]   .  CIII-62) 

-co  v        n  n 
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Thence 


iv  X 

*W  -  h(2K^.1A)  r  (1  -v;  -  )  ,        (in-63) 


and 


in  .X 


Now,  since 


s.(*)  =  /  s  (n)  dn  ,  (111-65) 

1       0 


then 


6  +00 

S,U)    =  ^-     /    [      /       77  tan  ^  elvnndv]    dn 
0        -00 


^-       /     i  tan  Mr  C    /      e     n      dn]   dv 
2lT   _co        V  2       0 


Therefore , 


S1(d>)    =   1(d))    .  (111-66) 


Similarly 


T1(c{,)    =   J(4»).  (IH-67) 


4H 


Thus , 

iv  X 


S1(X)  =  h(2Kn2-iA)  \    f  ^  -  e   n  3  •         ^II1-68) 
1     n-1   n 


and 


X 
Sp(x)  =  /  S  (n)  dn,  (111-69) 

0    x 


so 


X      iv  n 


1     n=l  n  0 


o      °°   -i       ,   iv  X 

Riray   =   f  {x  +  f [e   n   -  ™- 

i  n=l  n       n 

(111-70) 


etc.,  and  likewise  for  T„ 


C.   SOLUTION  THROUGH  COLLOCATION  TECHNIQUE 

Utilizing  the  same  geometry,  nondimensionalization,  and 
boundary  conditions  as  for  the  first  portion  of  Section  III, 
the  boundary  value  problem  is  solved  by  assuming  the  potential 
function  is  a  summation, 


*(X,Y)  =  $q(X,Y)  +  <j>0(X,Y)  +  <f>°(X,Y)  +  (})1(X,Y)  ,      (III-7D 


^5 


where  4>_  and  <$>     are  the  isolated  blade  potential  solutions 
and  <}>0  and  <f>  ,  are  the  interference  potentials  for  the 
reference  and  adjacent  blades. 

Such  an  approach  was  proposed  by  Gorelov  [Ref.  8]  for 
sub/supersonic  flows.   See  also  Samoylovich  [Ref.  21]. 

Utilizing  the  boundary  conditions  and  the  known  source 
distributions 


4>°y(X,0)  =  vm(X)  ,    m  =  0,1  ,  (IH-72) 


obtain  two  equations 


♦oy(X,Y0-0)  +  ^(X.Y^-p)  +  ^(X.Y^-p)  -  0 

(HI-73) 

Reference  Blade, 


*0Y(X'Y0=P)  +  WX>Y0=p)  +  *iy(X'Y1=0)  =  ° 

Adjacent  Blade.  (111-7^) 

The  isolated  blade  solution  is  already  known  from  Teipel 
[Ref.  30]  so 

K2        QY0 
n            i      X      Exp[%-(X-S)  -  ty^J 
*°(X,Yn)  -  -  -i-   /  vn(S)  S U"b;   dS,   (111-75) 

0       °         -yw  -i  -yx^ 


ne 


and 

.2 


2        QY 

0  1      X       Exp[^-(X-S)  -  -n^y] 

♦"(x.y.)  =  -A-   /  Vl(s)  S L2L±L_  dS.    (in-76) 

1  -1     "fQ1  -1    X  -yJZ^S 


It  is  expected  that  interference  potentials  will  have 
the  same  form  but  contain  unknown  source  distributions  un 

and  u, .      Thus , 

p 

2  Y 


,              X               Exp[^r-(X-S)    -    J  -^-^p   ]    dS 
♦  0(X,YQ)    =    -   J-        /   Uq(s)    9 J    (X-S) ^ 


-TO  -l  -Vx^s 


(111-77) 


and 

.2 


r2                         Y 
n              X               Exp[— (X-S)    -   J   (I      ,    dS 
<MX,Y..)    -  -i-       /  u,(s)    Ss IJXz^ 1         (111-78) 

1   -1   -^  -i  ^s" 


It  is  the  fundamental  point  of  this  solution  to  approxi- 
mate the  unknown  and  known  source  distributions  by  series, 
the  coefficients  of  which  are  to  be  determined  by  satisfying 
equations  (111-73  and  7*0  at  a  sufficient  number  of  points 
on  the  blades. 

Elder  [Ref.  7]  utilized  a  power  series, 


L-l 

u  (X)  =  Z  v  Xn  ,  (IH-79) 

m      n=Q  mn 
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and 


v(X)  =  I      e  xn  ,    m  =  0,1.  (111-80) 

m      n=0 


Elder  obtained  the  coefficients  of  these  series.   Expanding 
upon  his  methods,  this  author  obtained  potentials.   For  small 
numbers  of  collocation  points  (L<10),  A^O ,  K<1,  agreement 
with  the  potential  calculated  by  the  other  methods  was 
generally  good.   Numerical  difficulties  above  those  parameter 
ranges  and,  especially,  for  X=0  led  to  the  introduction  of 
the  ' step-wise-linear '  series  approximation 


L-l 
um(X)  =  vmn  +  vm1.X  +  I      vwir,(X-XM)H(X-Xr,)  ,         (111-81) 

n=2 


and 


L-l 
E 

n=2 


vm(X)  =  em0  +  6ml'X  +  L    9mn(X-Xn)H(X-Xn)  •         (III~82) 


where 


0  z  <  0 

H(z)  =  .  (III-83) 

1  z  >  0  . 


The  first  two  terms  were  necessary  to  allow  for  initial 
value  and  slope  of  the  source  distributions  (subsequent 
coefficients  represent  changes  in  that  slope  at  the  selected 
points,  XN).   Further,  they  allow  simple  representations  for 
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the  isolated  blade  source  distribution  coefficients.   For  a 
pitch  oscillation  about  the  midchord, 


Vx)  =  eoo  +  eoilX  > 


where 


600  =  -1  '     eoi  =  "  iK  »  (111-84) 


and 


v  <x)  =  e   +  e  .x  , 


where 


610  =  -  eia  ,    ell  =  -  iKeia.  (III-85) 


Inserting  these  step-wise-linear  approximations  into 
equations  (111-75,  76,  77,  and  78)  and  substituting  the 
resulting  expressions  into  the  boundary  condition  equations 
for  each  blade  (equations  III-73  and  7*0  obtain 
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? L  2 


V00   +    2*V?V10        ' 


^y^ — dS   Vqi  x 


2  2 

K    ,„   .v    Q   9 


(X-S)j/d 


-1 


n=yvon(x-vH(x-v 


X    (S-X_)H(S-X_)  K2 


+  2  #vm  _/    "3/2    ExpcVx-s)-f  P^pgj!  ds) 


--fff^io^ 


rX  Exp[^(X-S)-Q   4x3BT^ 


(X-S) 


372 


dS 


2  2 

+  e  /    s  S \  ^-^     ds] 

00    -1  (X-S)J/<d 


Reference  Blade,      (111-86) 

„       +  £  JffV  *  Exp[^(X-S)-a  ^   dS 

*10  +  I  Ji  Vo°  ./     ^372 +     u-   x 


+  f-/Fvoi./s 


?  2 

■K    ,„   0x    Q      p^ 


Exp[— (X-S)-f  T^y]    dS 

S A/ +      z    {Vl     (X-X    )H(X-Xj 

(X-S)J/^ 


n=2 


In  n  n' 


X    (S-X)H(S-Xn)  K2  Q      2 

+  £.  JS  v_        /     "      ^  n     Exp[K    (x_s)    Q   P-jdS} 


V°N  -/  (X-S)3/2 


•  K2 


Q  x  ¥    (X-S) 

2 


=    -    £  nE 


X  Exp[^-(X-S)-|   (P_S)]    dS 


2    Vtt 


[e  / 

00    -1  (X-S) 


372" 


rX  s  ^P[^(x-s)-f  pLst]  dS 
601  -i'    S  ^T72" 


(111-87) 


Adjacent   Blade. 
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These  equations  are  satisfied  at  each  collocation  point 
to  obtain  the  unknown  coefficients.   The  coefficients  can 
then  be  used  in  the  following  equation  for  the  potential, 

K2 

*(X'V°>  "  "  ^  "eoo  +W  ./        (x,Qs)i/2      dS 

K2 

X        Exp[~(X-S)] 

♦t.01*v01]^/  s_J--T75-dS 


X'(S-X    )H(S-X    )  K2 

+  n=2V°n  -1  (X-S)1/2   "     ExpCVX-S)]   dSJ 

t  2  2 

X  Exp[4r(X-S)-§.  ^r^y] 

+  —  ([enn+vin]   /  § ,  ,:  (X"S)  dS 

-iW   10  10  -l      (x-s)172^ 

2         2 
X   Exp[4-(X-S)-%  ?Y  oi] 

+  [e-.+v..]   /  s  § n/;  (X"S)  ds 

11  11  -1         (x-s)1/V 


X  (S-X  )H(S-X  )      v2  2 

+   E   V,     / S ~2-  Exp[%-(X-S)-H  P_.]dS} 

n=2   ln  -1      (X-S)1/2        Q       ¥  (X~S> 

(111-88) 


D.   EVALUATION  OF  THE  INTEGRALS  OCCURRING  IN  EQUATIONS 
(111-24,86,87  AND  88) 

The  integrals  encountered  are  of  the  form 


Y  XT  O  O 

In  =  -/    T^7"2"  ExpCVx-s)  -  rfcsy]  dS>        (III"89) 

N   =    0,1;      m  =   1,3. 
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Introducing  the  substitution,  U  =  X-S ,  separating  I 
into  its  real  and  imaginary  parts,  and  simplifying  as  did 
Elder,  obtain 


X+l  CX_U}N 
JNR  =   /      m/2   ExP^AU~^  cos[2KAU+^]  dU,         (111-90) 


and 


X+l  ,    ,N 

JNI  =  f  m/2   ExP^AU-g-]  sln[2KAU+|]  dU,        (III-9D 


where 


h   =  ^R  -  ±TNI    '    N  =  °^>  (IH-92) 

A  =  f   ,  ■   B  "  T  >   E  "  T  <IIJ-93) 


Then,  it  is  necessary  to  evaluate  six  integrals 

X+l 

=      J 
n 
I 


I   R  =      /  Uj"5/2   ExpCAAU-S]    °?!![2KAU+|]    dU,  (111-9*0 


where   j    =    1,2,3. 

These  are  four  cases  of  treatment  of  these  integrals: 

1.   The  Steady  Case 

Setting  K  =  0  reduces  the  integrals  to  a  relatively 
simple  case  which  may  be  evaluated  through  the  use  of  the 
error  function  or  by  numerous  numerical  techniques,  one  of 
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which  is  utilized  in  the  routines  described  in  the 
Appendices. 

2.   The  Oscillatory  Case  with  X  =  0 

In  this  case  the  integrals  are  of  the  form, 

Ic  -  /    UJ"5/2  sin[2KAU+|]  dU>  UII-95) 

s 

requiring  a  distasteful  integration  near  the  lower  limit  for 
J  <  5/2,  due  to  the  rapid  oscillations  and  increasing  magni- 
tude of  the  integrand  there. 

To  avoid  difficulties  at  the  lower  limit,  the  range 
of  integration  is  dissected,  i.e., 

X+l  XA 

j-5/2    coS|-2KAU+E-j    dU  =      f     ■  yJ-5/2    cos[2KAU+E-j      y 

0  sinL     u       0  sinL     u 


X+l 

+   /    UJ'~5/2  C?S[2KAU+|]  dU.     (111-96) 
XA  sinL      U 


The  limit  XA  will  be  chosen  so  as  to  allow  the  second 
integral  to  be  evaluated  numerically  while  not  greatly  af- 
fecting accuracy.   The  first  integral,  hereafter  the  'correc- 
tion integral',  will  be  approximated  in  the  manner  that 
follows . 

Recall 


cos[2KAU+-]  =  cos[2KAU]  cos[-]  -  sin[2KAU]  sin[]|]  , 
ij  u  u 
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and 


sin[2KAU+|]  =  sin[2KAU]  cosf|]  +  cos[2KAU]  sin[|]  , 


and,  the  truncated  series  expansion 

2 


cos[2KAU]  -  1  -  (2*AU) 


and 


(2KAU)3 


sin[2KAU]  -  2KAU  -  -^^ 


Then,  for  2KAU  sufficiently  small 


/    UJ'"5/2  C?S[2KAU+-]dU  =   /  '  UJ'"5/2  c?s[*±]  dU 
0  smL     u      0  smLu 

YA  3    YA 

*  (2ka)  /  uJ"3/2  slnr|]du  -  (2*A)   /  uJ""1/2  =?srf]du 

.~  cos   u  n  sin  u 


±   12K/03     r»  uJ+1/2   Si„  .     du_  (III.97) 

D      ,.  C  OS  U 


Since  integrals  of  the  form 


XA 

/   um/2  c?s[-3  dU  (111-98) 

.         sinLu 


5^ 


may  be  expressed  in  terms  of  the  Fresnel  Integrals,  the  Inte- 
grals in  (111-97)  can  be  calculated  after  integrating  by 
parts 

XA 


t        nJ-5/2    CQSpE.         =       2     xv3-3/2    cosr_E_n 
f        U  sinLUJdU        2J-3  sinLXAJ 


XA  U^3/2        sinrE 


+    2E      /       2 p   *„,[£]    dU>  (IH-99) 

0  (2J-3)TT   cos   L 


where,    j    >   3/2.      From  Elder    [Ref.    7]    for  J    =   1 


XA  rr^ 

/      u-3/2  °°*[§]dU  -  #  [|  -  cs(^)]  ,  (in-ioo) 


where 

E 
-                 XA    ^?S(t)    dt 
C(JL)    =      /        £in _  (III-101) 

s   XA  0  -V2lre 


Then,    for   j    =   2 

XA      _  /n  _  -  /0  -,  XATT-3/2    sinrE-,  .„ 

„-l/2    cosrEn,,,        „    ,r.-l/2  Jn  _    ^.„     ,  c0SL-JdU. 

(III-102) 


/  lT1/2  c?s[^]dU  =  2-XA"1/2  ^S[|rl  T  2-E  /  U     cosLU 
sin  u  sln  XA       » 


Thus,  the  correction  integrals  can  be  evaluated  sequentially 
on  increasing  j . 
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3.   The  General  Case,  K^O  and  A^Q 

In  this  case,  as  the  lower  limit  is  approached,  the 
exponential  function  in  the  integrand  is  depended  upon  to 
drive  the  integrand  magnitude  to  zero  more  rapidly  than  it 
is  magnified  by  the  negative  power  of  U.   Thus,  that  portion 
of  the  integral  near  zero  is  neglected  by  specifying  a  non- 
zero lower  limit  for  which  a  negligibly  small  error  is  pro- 
duced while  still  allowing  numerical  integration  to  be 
utilized. 

H.      The  Case  for  P=0,  j  =  2,  K^O  and  A^O 

For  this  case,  another  approach  for  that  portion  of 
the  integral  near  zero  must  be  taken.   Due  to  the  absence 
of  the  -B/U  in  the  exponential  function,  that  portion  cannot 
be  neglected  as  in  part  3.5  above. 

The  range  of  integration  is  again  dissected,  a 
'correction  integral'  is  defined,  and  XA  is  selected  simi- 
larly to  part  2.,  above.   The  'correction  integral'  integrand 
is  expanded  using  the  series  expansions  for  the  functions 
within  and  higher  order  terms  in  XA  neglected.   The  resulting 
series  is  integrated  term  by  term  and  programmed  for  correcting 
the  numerically  calculated  portion  of  the  integral  (see 
Subroutine  ZSUM  in  the  Appendices). 

E.   APPLICATION  AND  COMPARISON  OF  SOLUTIONS  TO  SEVERAL 
CASCADE  PROBLEMS 

In  order  to  gain  some  confidence  in  the  Collocation  method 

solution  before  extending  it  to  the  porous  walled  sonic  wind 
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tunnel  problem,  the  solutions  obtained  for  sonic  cascade 
interference  in  the  previous  sections  were  programmed  for 
digital  computation  (see  Appendices  A-C).   They  were  applied 
to  cascade  problems  differing  in  the  parameters  of  interblade 
phase  angle,  reduced  frequency,  Oswatitsch  parameter,  and 
interblade  distance  with  pitch  oscillations  about  the  cascade 
blade  midchord  position.   Many  of  those  results  are  shown  in 
this  section  to  illustrate  the  ranges  of  parameter  values  for 
which  agreement  between  solutions  can  or  cannot  be  expected. 
Because  it  is  the  object  of  this  section  to  illustrate  the 
basic  traits  of  the  solutions,  especially  those  of  the  Collo- 
cation method  solution,  agreement  will  be  considered  to  be 
a  5$  or  so  difference. 

The  complex  velocity  potentials  calculated  by  the  three 
methods  are  presented  in  three  sections,   First,  reduced 
frequency,  Oswatitsch  parameter  and  interblade  distance  were 
fixed,  the  interblade  phase  angle  varied.   Second,  interblade 
distance  and  phase  angle  were  fixed  and  the  remaining  parameters 
varied  to  illustrate  their  interrelationship  to  the  agreement 
between  the  Collocation  and  other  methods.   Third,  the  inter- 
blade phase  angle  was  fixed,  the  Oswatitsch  parameter  set  to 
zero  (the  governing  equation  is  thereby  linearized) ,  and  the 
reduced  frequency  and  interblade  distance  varied  to  illustrate 
the  difficulties  the  Collocation  method  had  satisfying  the 
linearized  transonic  potential  equation. 
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The  results  shown  indicate  that  the  Collocation  method, 
as  applied  here,  agrees  reasonably  well  for  lower  frequencies 
and  larger  Oswatitsch  parameter  with  the  other  methods.   For 
larger  frequency  and/or  small  Oswatitsch  parameter,  the  step- 
wise-linear  series  approximations  to  the  then  oscillatory 
source  distributions  become  oscillations  whose  amplitude  may 
grow  substantially  from  the  blade  leading  edge  aft.   For 
cases  in  which  the  amplitude  growth  is  large,  agreement  be- 
tween the  potential  calculations  of  the  Collocation  method  and 
the  others  deteriorates.   Furthermore,  increasing  the  number 
of  Collocation  points  fails  to  improve  agreement.   Finally, 
it  should  be  noted  that  when  agreement  was  poor,  the  source 
distribution  series  approximation  yielded  an  oscillatory 
result  whose  amplitude  growth  was  large. 

It  has  been  observed  that  a  step-wise-linear  series 
approximation  can  exhibit  this  characteristic  when  numerical 
errors  of  sufficient  magnitude  are  present.   Further,  the 
effect  of  increasing  the  number  of  Collocation  points  (and 
coefficients)  is  to  increase  the  growth  once  it  was  present, 
as  observed  above.   Numerical  errors  of  substantial  magnitude 
can  be  traced  to  the  calculation  of  the  integrals  (their 
values  are  highly  damped  oscillations  along  the  chord)  in 
the  Collocation  method  solution  above. 

The  failure  to  produce  accurate  potentials  when  an  oscil- 
latory series  approximation  was  obtained  can  be  attributed 
to  the  large  magnitude  of  the  coefficients  near  the  trailing 
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edge.   The  resulting  potential  calculations  will  then  involve 
the  difference  of  large  numbers. 

1.  Interblade  Phase  Angle  Variation 

Table  I  shows  real  and  imaginary  potentials  for 
K  =  .10,  A  =  .10,  p  =  1.0  and  a  =  tt,  tt/2,  and  0  radians 
(recall,  a  =  tt  corresponds  to  an  out-of-phase  oscillation  of 
the  blades)  for  the  Collocation,  Ftourier  Transform,  and  La- 
place Transform  solutions.   Further,  the  Collocation  and 
Fourier  methods  are  shown  for  K  =  .10,  A  =  .50,  p  =  1.0, 
and  c  =  0. 

Agreement  is  seen  to  be  fairly  good  for  a  >  0.   At 
a  =  0,  the  Collocation  method  fails  to  agree  with  the  others 
aft  of  about  X  =  -.70,  unable  to  fit  small  amplitude  imaginary 
part's  oscillation.   An  increase  in  thickness  (increased 
Oswatitsch  parameter)  improves  the  agreement  slightly. 

2 .  Frequency  and  Thickness  Variation 

Figure  3  shows  the  results  of  potential  calculations 
versus  position  along  the  cascade  blades  described  in  Section 
III-A  (-1  <_  X  <_  1),  for  Collocation  (step-wise-linear  series 
approximation),  Fourier  Transform,  and  two  cases  each  of  the 
Laplace  Transform  and  Collocation  (power  series)  methods. 
Interblade  distance,  rho,  was  fixed  at  one  and  interblade 
phase  angle,  a,  at  pi  radians  (out-of-phase  oscillation 
produces  greatest  interference).   Reduced  frequency,  K,  was 
given  values  .01,  .1,  .3,  .5,  1.0,  and  1.5.   The  Oswatitsch 
parameter,  A,  took  values  .05,  .1,  .2,  and  .5.   For  all  but 
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TABLE  I 

COMPARISON  OF  POTENTIAL  (REAL/IMAGINARY)  CALCULATIONS 
OF  COLLOCATION,  FOURIER,  AND  LAPLACE  SOLUTIONS  FOR 

VARIOUS  PHASE  ANGLES 


FOURIER  LAPLACE 


X 

COLLOCATION 

K  =  0.10, 

T  =  0.10,  p  = 

For  a=3 

.10  rad 

-.9 

0.391/  -0.946 

-.7 

1.048/  -2.571* 

-.5 

1.741/  -4.286 

-.3 

2.426/  -5.937 

-.1 

3.163/  -7.607 

0.1 

3.900/  -9.266 

0.3 

4.672/-10.920 

0.5 

5.459/-12.570 

0.7 

6.271/-14.220 

0.9 

7.104/-15.870 

For  a  =  1 

.60  rad 

-.9 

0.007/  -0.448 

-.7 

-0.497/  -1.574 

-.5 

-1.005/  -2.783 

-.3 

-1.503/  -3.927 

-.1 

-1.960/  -5.156 

0.1 

-2.446/  -6.299 

0.3 

-2.865/  -7.589 

0.5 

-3.340/  -8.670 

0.7 

-3.715/-10.110 

0.9 

-4.187/-11.000 

For  a  =  0 

,0  rad 

-.9 

0.505/  -0.064 

-.7 

0.503/  -0.029 

-.5 

0.498/  -0.037 

-.3 

0.507/   0.002 

-.1 

0.490/  -0.033 

0.1 

0.521/   0.047 

0.3 

0.466/  -0.052 

0.5 

0.563/   0.129 

0.7 

0.393/  -0.128 

0.9 

0.686/  0.296 

0.486/  -0.854  0.489/  -0.854 

1.146/  -2.525  1.150/  -2.527 

1.827/  -4.196  1.836/  -4.193 

2.530/  -5.860  2.540/  -5.861 

3.254/  -7.522  3.264/  -7.508 

4.000/  -9.180  4.017/  -9.176 

4.768/-10.834  4.782/-10.810 

5.558/-12.486  5. 580/-12 . 480 

6.370/-14.134  6.392/-14.110 

7.204/-15.780  7.209/-15.740 


0.092/-   0.404  0.094/   -0.404 

.0.425/   -1.605  -0.423/   -1.604 

•0.923/   -2.776  -0.918/   -2.776 

•1.410/   -3.956  -1.403/   -3.956 

.1.883/  -5.144  -1.872/   -5.142 

■2.344/  -6.342  -2.334/   -6.343 

.2.794/  -7.548  -2.776/   -7.544 

■3.228/  -8.764  -3.214/  -8.743 

■3.652/  -9.990  -3.631/  -9.987 

■4.064/-11.226  -4.0  39/-H.210 


0.542/  -0.009  0.545/  -0.009 

0.498/  -0.034  0.500/  -0.035 

0.499/  -0.025  0.501/  -0.025 

0.499/  -0.015  0.501/  -0.01.6 

0.499/  -0.005  0.501/  -0.006 

0.499/   0.005  0.501/   0.005 

0.499/   0.015  0.501/   0.014 

0.499/   0.024  0.501/   0.024 

0.499/   0.0?4  0.502/   0.034 

0.499/   0.044  0.503/   0.045 
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TABLE   I    (Continued) 


X 

COLLOCATION 

FOURIER 

LAPLACE 

K   =    0.10, 

T    =    0.50,    p    =    1 

,00 

For  a   =   0, 

0   rad 

-.9 

0.43V  -0.088 

0.437/   -0.088 

-.7 

0.509/  -0.042 

0.500/   -0.039 

-.5 

0.489/  -0.024 

0.499/   -0.027 

-.3 

0.515/  -0.023 

0.499/   -0.017 

-.1 

0.485/      0.001 

0.499/    -0.007 

0.1 

0.520/   -0.009 

0.499/      0.003 

0.3 

0.480/     0.030 

0.499/      0.013 

0.5 

0.525/      0.000 

0.499/      0.023 

0.7 

0.476/      0.063 

0.499/      0.033 

0.9 

0.528/      0.003 

0.499/      0.043 
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Fig.  3-c,  in  which  twenty  were  used,  ten  potentials  were 
calculated  along  the  blade  and  those  positions  were  used  as 
node  points  for  the  Collocation  routines. 

Fourier  results  were  machine  plotted  as  a  solid  curve 
connecting  individual  points,  which  can  result  in  a  jagged 
curve  for  some  parameter  combinations  because  of  plotter 
limitations.   Collocation  (step-wise-linear)  results  are 
plotted  as  crosses.   Laplace  Transform  calculations  were 
made  for  K  =  .01,  A  =  .05  (computer  time  usage  =  6  minutes) 
and  K  =  1.0,  A  =  .1  (usage  =  HO   minutes)  and  plotted  as 
squares.   Power  series  Collocation  results  were  also  plotted 
as  squares  for  K  =  A  =  .1  (usage  =  1^  minutes)  and  K  =  .5, 
A  =  .1  (usage  =  50  minutes).   Typical  times  for  similar  cal- 
culations by  Fourier  routine  are  2-5  seconds  and  step-wise- 
linear  Collocation,  10-20  seconds.   Interestingly,  though  not 
surprising  considering  the  method  by  which  integrals  were 
calculated,  the  Collocation  results  for  20  points  were  obtained 
in  slightly  less  time  than  10  points. 

Figure  3  is  arranged  in  order  of  ascending  frequency 
with  A  varying  at  each  to  show  the  characteristics  of  the  Collo- 
cation method  solution.   Generally  good  agreement  in  the  cal- 
culation of  the  interference  potentials  for  all  four  methods 
used  is  shown.   Because  Fourier  and  Laplace  agree  very  well 
at  the  points  shown  and  for  many  others  not  shown,  the  Fourier 
solution  was  concluded  to  be  a  good  standard.   It  was  further 
concluded  that  increasing  frequency  adversely  affects  the 
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accuracy  of  Collocation  results  and  increasing  thickness,  as 
represented  by  the  Oswatitsch  parameter,  improves  it.   Further- 
more, increasing  Collocation  points  does  not  necessarily 
improve  accuracy  as  illustrated  by  Fig.  3-c  compared  to 
Fig.  3-b. 

It  was  observed  that  in  those  cases  for  which  rela- 
tively poor  potential  agreement  between  the  Collocation  and 
the  other  methods  was  obtained,  the  step-wise-linear  series 
source  distribution  approximation  was  an  oscillation  whose 
amplitude  grew  toward  the  trailing  edge.   The  coefficients 
of  those  series,  because  they  represent  the  changes  in  the 
slope  of  the  approximation,  can  illustrate  the  point.   For 
example,  for  K  =  .3,  A  =  .05,  which  shows  good  agreement,  the 
first  and  last  coefficients  are  (real/imaginary)  -.7280/1. 760 
and  1. 488/-. 0269,  respectively.   For  K  =  .5,  X   =  .05,  poorer 
agreement,  with  -.1357/1.879  and  -1.66V-2.942.   For  K  =  1.0, 
X  =  .1,  very  poor  agreement,  with  3-549/1.402  and  -164. V-392 . 7 
For  K  =  .3  and  A  =  .1  and  20  collocation  points,  poor  agree- 
ment, with  -.4969/3.505  and  1714 .0/-2466 .0.   The  same  result 
is  seen  in  the  reduction  in  accuracy  shown  for  a  =  0.0  in 
the  previous  section. 
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Fig.  3-a.   Complex  Velocity  Potential  Comparisons  for  o--n 
and  p=1.0  with  Varying  Frequency  and  Oswatitsch  Parameter. 
Velocity  Potential  vs.  X  (1  unit/inch). 
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Fig.  3b.   Complex  Velocity  Potential  Comparisons  for  o  =  tt 
and  p=1.0  with  Varying  Frequency  and  Oswatitsch  Parameter. 
Velocity  Potential  vs.  X  (1  unit/inch). 
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Pig.  3c.   Complex  Velocity  Potential  Comparisons  for  c  =  tt 
and  p=1.0  with  Varying  Frequency  and  Oswatitsch  Parameter. 
Velocity  Potential  vs.  X  (1  unit/inch). 
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Fig.  3d.   Complex  Velocity  Potential  Comparisons  for  a  =  ir 
and  p=1.0  with  Varying  Frequency  and  Oswatitsch  Parameter. 
Velocity  Potential  vs.  X  (1  unit/inch). 
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Fig.  3e.   Complex  Velocity  Potential  Comparisons  for  a=ir 
and  p=1.0  with  Varying  Frequency  and  Oswatitsch  Parameter. 
Velocity  Potential  vs.  X  (1  unit/inch). 
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Fig.  3f.   Complex  Velocity  Potential  Comparisons  for  a=ir 
and  p=1.0  with  Varying  Frequency  and  Oswatitsch  Parameter. 
Velocity  Potential  vs.  X  (1  unit/inch). 
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3.   Ostoatitsch  Parameter  =0.0 

Fig.  4  is  a  compilation  of  selected  interference 
potential  calculations  for  X   -   0.0.   Because  the  accuracy  of 
the  correction  integral  technique  used  in  the  Collocation 
method  (and  Laplace  Transform  method)  (see  Appendix  A-II)  , 
connects  accuracy  to  frequency  and  interblade  spacing,  various 
interblade  distances  were  chosen  to  illustrate  the  effect  of 
frequency  on  the  Collocation  method  solution. 

Again,  the  values  of  potential  are  plotted  against 
position,  X,  along  the  cascade  blade  (-1  <_  X  <   1)  ,  Fourier 
Transform  results  are  plotted  as  a  solid  line  connecting 
points,  Collocation  (step-wise-linear  series)  results  are 
plotted  as  crosses  (except  at  p  =  12,  where  Teipel  results 
for  an  isolated  wing  are  so  plotted) ,  and  Laplace  Transform 
results  are  plotted  as  squares.   The  number  of  collocation 
points  varied  from  ten  to  twenty,  as  can  be  seen  from  the 
plots. 

Fig.  4-a  shows  results  for  interblade  distance  rho 
equal  to  one.   For  K  =  .01,  the  Collocation  results  were  so 
bad  the  real  parts  would  not  plot  on  the  graph.   For  K  -  .1, 
the  real  parts  of  the  potentials  agreed  fairly  well.   The 
bottom  plots  show  the  effect  of  increasing  the  number  of 
points.   It  is  seen  again  that  the  Laplace  and  Fourier 
Transform  solutions  agree  rather  well,  as  indicated  in  the 
few  results  shown  and  others  not  shown. 
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Fig.  J|-b  and  c  show  results  for  p  =  H    (p  is  increased 
so  that  frequency  can  be  increased  while  still  attempting  to 
maintain  accuracy  in  the  Collocation  routine,  see  Section 
III-D  and  Appendix  A-II).   Agreement  is  only  fair  and  the 
conclusion  must  be  drawn  that  the  derivatives  of  potentials 
needed  for  calculations  of  pressure  coefficients,  lifts  and 
moments  are  likely  to  be  very  questionable.   As  in  the  pre- 
vious part  of  this  section,  large  growth  of  the  source  dis- 
tribution approximation  amplitude  would  signal  poorer 
agreement  for  the  Collocation  method. 

The  bottom  graphs  of  Pig.  4-d  are  for  p  =  12.0.   The 
crosses  depict  potentials  calculated  by  a  routine  following 
the  work  of  Teipel  for  an  isolated  wing.   The  convergence  of 
Fourier  Transform  solution  results  into  isolated  wing  values 
for  increasing  interblade  distance  will  be  discussed  later. 
Here,  these  data  further  serve  to  illustrate  the  applicability 
of  the  Fourier  Transformation  solution. 
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Fig.  4a.   Complex  Velocity  Potential  Comparisons  for  c=it 
and  X=0.  Varying  Frequency  and  Interblade  Distance. 
Velocity  Potential  vs.  X  (1  unit/inch). 
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Fig.  4b.   Complex  Velocity  Potential  Comparisons  for  o  =  tt 
and  X=0,  Varying  Frequency  and  Interblade  Distance. 
Velocity  Potential  vs.  X  (1  unit/inch). 
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Pig.  4c.   Complex  Velocity  Potential 'Comparisons  for  a=ir 
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IV.   TRANSONIC  CASCADE  AND  BIPLANE  STABILITY 

A.   GENERAL 

In  the  previous  sections,  the  Fourier  Transform  solution 
has  emerged  as  a  standard.   The  method  has  been  made  flexible 
by  the  inclusion  of  the  Oswatitsch  parameter  and  variable 
interblade  phase  angle.   It  has  been  compared  to  Laplace 
Transform,  Collocation  method,  and  Teipel  [Ref,  30]  solutions 
for  many  different  sets  of  parameters  and  has  produced  good 
results. 

The  Fourier  Transform  solution  has  been  applied  to  the 
problem  of  transonic,  unstaggered  cascade  aerodynamic  stabili- 
ty.  Selected  results  are  included  to  show  the  effects  of  fre- 
quency, blade  thickness,  interblade  distance,  and  interblade 
phase  angle. 

The  biplane  problem  was  solved  by  utilizing  the  Fourier 
solution  for  the  interference  caused  moment  upon  the  upper 
surface  of  the  lower  blade  and  adding  the  Teipel  [Ref.  30] 
solution  for  the  isolated  wing  moment  on  the  lower  surface 
of  that  blade. 

Isolated  wing  lift  and  moment  results  published  by  Nelson 
and  Berman  [Ref.  17],  Teipel  [Ref.  30],  and  Hamamoto  [Ref.  10] 
were  utilized  for  comparisons.   Barton  [Ref.  2]  and  Rott 
[Ref.  21]  were  also  checked.   The  Fourier  Transform  solution 
reproduced  all  of  the  aforementioned  results  with  reasonable 
accuracy,  even  when  the  interblade  distance  was  about  20 
(necessary  for  comparisons  when  A  =  0.0). 
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3.   TRANSONIC,  UNSTAGGERED  CASCADE  STABILITY 

Calculations  of  stability  boundaries  were  made  for  certain 
interblade  phase  angle,  a,  blade  thickness/chord  ratio,  x, 
and  interblade  distance,  p.   The  moment  coefficients  were 
then  calculated  for  varying  frequency  to  find  the  frequency 
at  which  the  moment's  phase  angle  was  -  180  or  0,  for  each 
pivot  location.   Blade  thickness  was  converted  to  Oswatitsch 
parameter  as  per  Teipel  [Ref.  30], 


A  =  1.559  [(y+Dt]2/3  (IV-1) 


The  boundaries  are  plotted  so  that  stable  ranges  lie  above 
the  curves  which  represent  the  boundaries  of  stability  for 
that  p. 

It  is  seen  that  for  p  =  1,  the  boundary  frequency  lies 
off  the  curves  for  forward  pivot  locations.  Further,  fre- 
quency, interblade  distance,  and  blade  thickness  increases 
were  stabilizing  as  were  interblade  phase  angle  reductions 
and  rearward  pivot  shifts.  It  should  be  noted  that  the  effect 
of  interference  was  destabilizing  as  indicated  by  the  strong 
effect  of  decreasing  interblade  distance,  p. 
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C.   TRANSONIC  BIPLANE  STABILITY 

A  single  result  (Pig.  5-f)  was  included  to  show  the  effect 
of  interference  on  the  Biplane.   Because  the  blades  have  Iso- 
lated blade  pressures  on  their  'outer'  surfaces,  the  effect 
of  interference  is  much  less  evident  than  observed  in  cascade 
flows.   This  can  be  seen  by  comparing  the  cascade  and  biplane 
curves  for  T  =  .01,  and  a   =  7r(Figs.  5-c  and  5-f).   For  close 
interblade  spacing,  where  interference  is  most  prevalent,  the 
biplane  stability  curves  drop  to  lower  frequencies  at  more 
forward  pivot  locations,  indicating  greater  stability. 
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V.   EXTENSION  OF  THE  COLLOCATION  SOLUTION  TO  THE 
POROUS  WALLED  TRANSONIC  WIND  TUNNEL 

A.   PROBLEM  FORMULATION 

1.   Geometry  and  Boundary  Conditions 

The  geometry  considered  here  will  differ  in  that  the 
adjacent  blade  of  the  Collocation  method  solution  for  the 
sonic  cascade  is  now  a  flat  wall,  the  porosity  of  which  may 
be  varied. 

The  boundary  condition  acting  at  the  reference  blade 
(or,  now,  the  airfoil)  is  as  before 


8h 
if,v(X,Yn=0)  =  vn(X)  =  -j^  +  iKhu.  (111-4) 


From  Drake    [Ref.    6],    the    condition   at   the   wall    can 
be   stated  as 


i|>Y(X,Y=p)    +    awj![iKiKX,Y  =  p)    +    ij;x(X,Y=p)]    =    0.  (V-l) 

The  condition  of  zero  upstream  disturbance  holds  as 
previously.   Furthermore,  non-dimensionalizations  are  the 
same . 

2.   Solution  of  the  Oscillating  Sonic  Flow  Equation 

As  in  Section  III,  the  source  distributions  shall  be 
assumed  to  be  series.   However,  because  the  'adjacent  blade' 
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is  now  a  flat  wall,  Its  'isolated  blade'  potential  and 
derivatives  vanish.   Thus, 


*(X,Y0)  =  *o(X,Y0)  +  *0(X,Y0)  +  cf,1(X,Y1).  (V-2) 


Now,  from  the  boundary  conditions  and 


<J>°Y(X,Y0=0)  =  vm,    m  =  0,1;  (V-3) 


the   airfoil   and  wall   equations   become 


♦oy(x>V0)  +  4>1y(x>y:l=-p)  =  °»  (v_i,) 


and 


<fr°Y(X,Y0  =  p)  +  (t>0Y(x,Y0  =  P)  +  ^(x.y^o)  + 


<W[C(X>Vp)    +    WX>Vp)    +    *ix<X>V0)]   + 


iK[<J>°(X,Y0=p)    +    ({,0(X,Y0  =  p)    +    (|)1(X,Y1=0)]}    =    0 


(V-5) 


3.   Derivation  of  Expressions  for  Derivatives  of  the 
Potential 

In  Section  III-C  the  isolated  and  interference  poten- 
tials are  expressed  by  integrals  of  the  appropriate  source 
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distributions  (see  equations  111-71,  75,  76,  77,  and  78). 
The  wall  boundary  condition  requires  derivatives  on  the  X 
variable  of  those  potentials.   Thus,  the  following  derivatives 
must  be  evaluated, 

DIN  =  »  '   /  S»  E^JiX-S)    -  S  ^  «.       (V-6) 
dA  XN  l/x^ 

n  =  0,1;   P  >  0. 


For  n  =  0, 


2  2 

.       X  EXP[^-(X-S)  -  S   -TY-Try] 

DI0  =  3    / 8 2 U"b}   dS.  (V-7) 

dA  XN  VX^S1 


For  n  =  1, 

2  2 

DI1  =  ^(X    /   9 4  U"S;  dS 

dA    XN  i/X^S1 

2  2 

_   /  9 ?/p     dS}  •    <v-8) 

XN  (X-S)  U 


Thus,  derivatives  are  required  of  two  integrals 


DIM  =  -^        /    (X-S)mExp[^-(X-S)    -   \  xfcsT]dSj      m  =   ±   \ 
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Applying  Leibniz'  rule 
£  /qf(s,a)dS  =  /  £f(s,a)ds  +  f(q.a)g  -  f(P,a)g  , 

DIM  becomes 


DIM  =  xn^  ^<(x-s^^(x-s)  ~  *  Tx4sr]}dS 


+  (x-s)mExP[^-(x-s)  -  f  T|rsy^l     • 


To   differentiate   the    integrand,    let   U  =    X-S, 


-^{UmExp[AU   +   B/U]}    = 


{mUm   -1  +   K2   TTm        QP2   TTm-2}   _      .K2  Q      P2-, 

_  U      +   -jp  U  Exp[-Q-U   -   ^  -j^] 


There   are  two   cases   to   consider,   m  =   ±1/2   that   is 

/->  X—X  /-,  ■^•— X„  r-    /O 

«-2  n        ,  /0  n   2  n  TT-5/2 

DI-   =   ~     /  U_1/^Exp[C]dU  +  =6-     /  U  exp[C]    dU 

y    0  M      0 

1        X~Xn     _o/2                            Exp[-S  L-]    0 
-  4     /  U    J/^Exp[C]    dU   +   TTTL-\       »         (V-10) 

^  o  ir'^ 


and 
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2  n  0 

DI+   =  Jr.      /  U"1/2Exp[C]dU   +   SL-     /  lT3/2Exp[C]dU 

Q    0  4      0 

n  2       0 

+   |      /  '   U~1/2Exp[C]dU   +   U1/2Exp[-^     \l\       ,         (V-ll) 


where 


2  2 

c  -  51  u  -  Q  p 


The  last  term  in  equation  (V-ll)  vanishes  nicely  for 
all  P.   As  long  as  P  >  0,  the  last  term  in  (V-10)  does  as 
well.   For  P  =  0,  however,  care  must  be  exercised. 

For  P   =    0, 

r\  A""  A.  <2 

DI-(P=»0)    =  5r-     /  U~1/2Exp[^4J]dU 

q  0  y 

-  |     /  U~3/2Exp[~U]dU   +    -^-1       . 

0  U 

* 

Since 


ax  ,        ax  ax 

'      m  m-1     m-1  m-1      J^-l        ' 

X  XX 


then 

2 
X~Xn        Q/5  K2  ExpC^U]   X"Xn 

/  U~3/2Exp[~U]dU   =   -2 Jfc-\ 

0  y  V  '        0 

2 

K2        X~Xn  Exp[%-U] 
+   2  |-     /  j4—  dU. 
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Thus, 


K2 
Exp[^-(X-X  )] 

DI-(P-O)  =  g     n    .  (V-12) 
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B.   RESULTS  OF  POROUS  WALL  TUNNEL  CALCULATION 

The  results  of  the  preceding  analytical  work  were  pro- 
grammed (see  Appendix-D)  and  pressure  coefficients  calculated 
Comparisons  were  made  with  solid  wall  tunnel  and  sonic  free 
jet  results  obtained  from  the  other  methods.   Experimental 
comparisons  were  made  in  the  steady  flow  case  using  results 
of  Spreiter,  Smith  and  Hyett  [Ref.  27]. 

Figures  6,7,  and  8  plot  solid  wall  tunnel  program  results 
as  crosses,  sonic  free  jet  as  squares,  and  porous  wall  pro- 
gram results  as  solid  lines  connecting  data  points.   Further, 
the  porosity  is  printed  to  the  right  of  each  curve  of  collo- 
cation method  results.   Results  from  [Ref.  27]  and  for  a  run 
of  the  Collocation  method  for  variable  porosity  are  plotted 
as  diamonds  and  triangles  on  Figs.  6  and  7,  respectively. 

The  results  from  [Ref.  27]  plotted  on  Fig.  6  are  for  a 
circular  arc  airfoil  in  a  choked  wind  tunnel  with  no  wall 
porosity.   The  agreement  is  relatively  good  but  it  must  be 
considered  that  the  use  of  a  constant  fluid  acceleration  (X) 
is  an  approximation.   Further,  the  experimental  results  indi- 
cate a  magnitude  reduction  in  the  pressure  coefficients  with 
introduction  of  porosity  in  the  manner  shown  by  the  results 
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of  the  computer  routine.   Porosities  used  were  0.0,  .05, 
.5,  10.0,  and  100.0. 

Figs.  7  and  8  are  the  results  of  nonsteady  analyses  for 
which  there  is  no  experimental  or  theoretical  comparison 
available.   It  can  be  mentioned,  however,  that  for  both  the 
increases  of  porosity  produces  the  expected  trend  in  the 
porous  wall  tunnel  collocation  routine  results.   Note  that 
the  pressure  coefficients  shown  are  normalized  to  the 
magnitude  of  the  pitch  oscillation,  aQ . 

It  should  be  noted  that  for  K  =  .5  (Fig.  8),  the  Collo- 
cation method  solution  was  near  the  limit  of  its  applicability, 
as  expected  from  the  results  of  Section  III-E.   In  particular, 
for  a  0  =  0,  the  source  distribution  series  approximation  was 

Wj6 

a  mildly  divergent  oscillation  and  the  result  shown  in  Fig. 
8-a  is  oscillatory  as  well.  Interestingly,  as  porosity  was 
increased  the  character  of  the  results  improved. 

The  Collocation  routine  was  modified  for  the  utilization 
of  the  following  variable  porosity, 


avl   =  5.0  (X  +  1.0) 


Thus,  the  porosity  varies  linearly  from  zero,  at  the  leading 
edge,  to  ten  at  the  trailing  edge.   Of  course,  there  is  pre- 
sently no  way  to  verify  the  results  but  they  are  promising. 
However,  the  source  distribution  series  approximation  was 
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more  oscillatory  for  the  variable  porosity  case  than  that  for 
constant  porosity,  suggesting  further  analysis  is  necessary 
to  confirm  the  Collocation  method's  applicability  to  this 
technique. 
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Pig.  6.  Steady  Flow  Pressure  Coefficient  vs.  X,  t  =  .06, 
for  Various  Porosity  (0. , .05,  .50,10.  ,100. )  and  Solid  Wall 
Tunnel,  Sonic  Free  Jet,  and  Experimental  Results. 
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Fig.  7a.   Nonsteady  Flow  Normalized  Pressure  Coefficient 
vs.  X,  K  =  .1,  t  =  .06,  for  Various  Porosity  (0.  ,  .05,  .50,5.0, 
100.),  Solid  Wall  Tunnel,  Sonic  Free  Jet,  and  Variable 
Porosity  Collocation  Routine  Results.   Real  Part. 
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Fig.  7b.   Nonsteady  Flow  Normalized  Pressure  Coefficient 
vs.  X,  K  =  .1,  t  =  .06,  for  Various  Porosity  (0.,.05,.50, 
5.0,100.),  Solid  Wall  Tunnel,  Sonic  Free  Jet,  and  Variable 
Porosity  Collocation  Routine  Results.   Imaginary  Part. 
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Fig.  8a.   Nonsteady  Plow  Normalized  Pressure  Coefficient 
vs.  X,  K  =  .5,  t  =  .06,  for  Various  Porosity  (0.,. 05, 1.0, 
5.0,100.),  Solid  Wall  Tunnel  and  Sonic  Free  Jet  Results. 
Real  Part. 
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Fig.  8b.   Nonsteady  Flow  Normalized  Pressure  Coefficient 
vs.  X,  K  =  .5,  t  =  .06),  for  Various  Porosity  (0., .05,1.0, 
5.0,100.),  Solid  Wall  Tunnel  and  Sonic  Free  Jet  Results. 
Imaginary  Part . 


97 


C.   DISCUSSION  OF  MURASAKI  AND  DRAKE  POROUS  WALL  SONIC  WIND 
TUNNEL  SOLUTIONS 

Murasaki  [Ref.  18]  solved  the  steady  sonic  flow  porous 
wall  wind  tunnel  problem.   His  results  could  not  provide  more 
than  one  check  case  against  the  Collocation  method,  though. 
However,  they  show  his  solution  to  have  good  agreement  with 
those  of  experiment  for  a  solid  wall  tunnel  in  [Ref.  27] 
(those  results  are  shown  in  Fig.  6) . 

Drake  [Ref.  6]  Laplace  transformed  the  nonsteady  linear- 
ized potential  equation  and  applied  the  porous  wall  wind 
tunnel  boundary  conditions,  leading  to  the  following  inversion 
for  the  Green's  function, 

(X  0+)  =  "T""1  ^B  cosh[^Ikp  p]  +  q(p+iK)  sinh[^2lKl?  p]     ( 
*-     jlfeiKp  [-^THp*  sinhftf2lKp  p] +a(p+iK)  coshEtfTiKi?  p]]J 

Drake's  results  of  lift  coefficient  could  not  be  used  as 
a  verification  because  of  the  known  inability  of  the  Collo- 
cation method  to  produce  derivatives  of  the  velocity  potential 
for  the  linearized  (zero  thickness  effect)  case.   Thus,  Drake's 
solution  was  modified  to  include  the  Oswatitsch  parameter  and 
computer  programmed.   The  results  of  that  program  did  not 
agree  with  reasonable  estimates  so  the  Drake  solution  method 

was  followed  to  produce  a  solution  for  a    „  =  0 ,  in  order  to 

wet- 
test the  method  at  as  simple  a  level  as  possible. 
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For  a   ,    =    0,    the   potential   is,    using  notation   of  Section 


III-B 


+  X  + 

<KX,0    )    =      /   a(t)    g(X-t,0    )    dt, 

0 


where 


a(t)    =    1   +    iK(t    -   A)    , 

g(X,0+)    =   |i.  fy+±C°  e^X   coshi^  P^        dp 

Y-i00  VQf?  sinhfTOp1  p] 

L-   •     .     s  x 

=    e^      E      4  e   n      , 
n=0    pQ 


and 


Sn  =  -H-5-  .,        n=   0,1,2,... 

Qp 


Though  this  result  is  derived  in  a  parallel  manner  to 

Drake's  a    .    ?   0,  A  =  0  result  (inversion  by  the  Theorem  of 

Residues),  it  leads  to  a  series  whose  convergence  is 

questionable,  i.e., 

2  2 
n  tt  v 

oo      -   — -  X 

I     e  QP2   . 
n=0 
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In  fact,  computer  generated  results  using  the  above  for- 
mulation showed  good  agreement  with  previous  results  of  Four- 
ier and  Laplace  Transform  and  Collocation  method  routines  for 
K  <  .3  and  X  ?   0.   In  those  cases  for  which  agreement  was  not 
obtained,  it  is  believed  that  the  computer  routine  signaled 
series  convergence  at  a  local  minimum  of  term  value. 

It  was  concluded  that  a  Laplace  Transform  solution  of 
this  type  needs  further  analysis  before  it  could  be  used  in 
a  verification  of  Collocation  results. 
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VI.   DISCUSSION  AND  SUMMARY 

A.   ANALYSIS  OF  SOLUTIONS 

Reference  should  be  made  to  Section  III-E. 

1.  Laplace  Transform 

This  solution  provided  a  reliable,  if  time  consuming, 
check  on  the  others.   Despite  the  principal  shortcoming  of 
an  unsophisticated  integration  technique,  it  always  gave 
reasonable  agreement  with  the  Fourier  Transform  solution. 
Even  for  X  =  0,  the  case  that  required  the  'correction  inte- 
gral' described  in  Section  III-D,  agreement  with  the  Fourier 
solution  was  good  for  the  parameters  tested.   Thus,  the  cal- 
culation of  the  integrals  utilized  in  the  Collocation  method 
solution  was  verified. 

As  programmed,  this  solution  was  too  extravagant  in 
its  computation  time  usage,  principally  due  to  the  time  re- 
quired for  integrating.   Presently,  integrals  are  calculated 
from  the  leading  edge  to  the  X  value  for  which  a  potential  is 
desired.   Perhaps,  a  better  approach  would  be  to  calculate 
the  integrals  in  increments  between  each  X  value,  as  does  the 
Collocation  routine  (see  Appendix-C).   Further,  an  improved 
integration  routine  itself  may  speed  the  results  and  improve 
their  accuracy. 

2.  Fourier  Transform 

The  Hamamotc  Fourier  Transform  method,  as  extended 
here  for  variable  blade  thickness  (via  the  Oswatitsch 
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parameter),  frequency  (to  Include  K  =  0) ,  and  interblade 
phase  angle,  was  a  reliable  and  reasonably  rapid  method  for 
velocity  potential,  pressure  coefficient,  lift,  and  moment 
calculations  (which  includes  the  solid  wall  tunnel  and  free 
jet  problems).   It  must  be  noted  that  for  large  interblade 
distances,  beyond  those  required  to  initially  obtain  isolated 
wing  results,  the  series  involved  may  not  converge.   This 
must  be  considered  principally  for  the  A  =  0  case. 
3.   Collocation  Method 

The  Collocation  method,  as  programmed  here,  was  a 
qualified  success  due  to  its  agreement  limitations  for  higher 
frequencies  and  lower  Oswatitsch  parameter,  as  discussed  in 
Section  III-E.   Whenever  the  other  methods  showed  that  the 
potential  results  were  not  reasonably  smooth  or  large  magni- 
tude differences  were  present  between  the  real  and  imaginary 
parts  of  the  potential,  the  coefficients  of  the  series  became 
large  and  the  computed  velocity  potentials  deviated  from  those 
of  the  Fourier  and  Laplace  Transform  solutions.   Increasing 
the  number  of  collocation  points  did  not,  in  general,  improve 
the  agreement.   It  is  believed  that  the  character  of  the  step- 
wise-linear  series  approximation  in  the  presence  of  the  numeri- 
cal inaccuracies  of  the  integral  calculations  (though  not 
terribly  significant  in  themselves)  can  lead  to  large  errors 
in  the  calculated  potential,  especially  near  the  trailing  edge 

Because  of  the  inability  of  the  method  to  reliably 
produce  potentials  for  A  =  0,  it  cannot  calculate  pressure 
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coefficients  for  that  case,  which  are  needed  if  lift  and 
moment  results  are  desired. 

Perhaps  an  improvement  to  the  step-wise-linear  series 
is  the  Fourier  series.   It  might  overcome  the  numerical  pit- 
falls encountered  by  the  power  series  approach  of  Elder  and 
the  aforementioned  qualification  of  present  series  by  pro- 
viding a  better  approximation  to  the  source  distribution  when 
more  complex  behavior  is  expected. 

B.   CASCADE  STABILITY 

Two  qualifications  must  be  mentioned  as  to  the  applica- 
bility of  these  results,  produced  by  the  Fourier  Transform 
solution.   First,  the  t  =  0  results  were  a  solution  to  the 
linearized  potential  equation,  the  applicability  of  which  is 
a  function  of  reduced  frequency,  K,  e.g.  Landahl  [Ref.  13] 
stated  the  following  requirement  for  applicability, 

K  >>  |1  -  M.  | 

where,  MT  =  local  mach  number. 
JLi 

Secondly,  for  t  ¥■   0,  the  approximation  of  the  Oswatitsch 
parameter,  X,  was  utilized.   This  can  only  be  a  rough  estimate 
to  the  real  behavior  but,  at  least,  it  is  an  initial  step. 
A  more  accurate  analysis  could  perhaps  be  made  through  the 
application  of  Spreiter's  local  linearization.   (See  part 
D  of  this  Section) . 
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C.  POROUS  WALL  WIND  TUNNEL  INTERFERENCE 

The  Collocation  method  applied  to  the  porous  walled  wind 
tunnel  is  subject  to  the  same  qualifications  discussed  earlier 
in  this  section  for  the  Collocation  method. 

The  results  are  difficult  to  verify  experimentally  due 
to  the  paucity  of  such  data.   Further,  theoretical  verifica- 
tion is  difficult  to  obtain  as  only  Murasaki's  approach  is 
available  for  steady  flow  and  Drake's  analysis  for  nonsteady 
porous  wall  tunnel  flows.   Unfortunately,  Murasaki's  data 
could  be  used  to  verify  only  a  single  steady  flow  case  and 
the  Drake  Laplace  Transform  solution  was  restricted  to  the 
linearized  case  for  which  the  Collocation  method  solution 
would  not  produce  reliable  results.   However,  the  trends  of 
the  Collocation  method  solution  are  correct  and  the  results 
for  very  large  and  very  small  porosities  compared  well  with 
the  free  jet  and  solid  wall  results. 

The  modification  of  the  routine  for  variation  of  porosity 
along  the  wall  appears  promising  and  may  allow  a  numerical 
test  of  the  Sears  [Ref.  25]  technique.   Further,  a  kind  of 
local  linearization  may  be  applied  to  collocation  methods 
(see  the  next  part  of  this  section). 

D.  FUTURE  WORK 

Various  areas  of  improvement  have  been  recognized. 
1.   Improve  series  approximation  methods  for  Collocation 
methods . 
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2.  Improved  numerical  integration  schemes  for  those  methods 

3.  Application  of  Spreiter's  'local  linearization'  to 
Fourier  Transform  and  Collocation  methods. 

The  first  two  have  been  discussed.  The  third  requires 
the  following  elaboration  and  description  for  steady  sonic 
flow. 

Recall  that  the  Oswatitsch  parameter  was  used  to  approxi- 
mate the  nonconstant  coefficient  of  the  governing  equation, 
i.e . , 


became , 

where,  A  =  constant  >  0,  assuming  that  the  fluid  acceleration 
was  slowly  varying. 

The  solution  for  the  steady  flow  potential  on  a  single 
foil  is, 

7(X,0)  =  —   /  — -^-  dS. 

■yo  -1  -Jx^S 

Now,  Spreiter  argued  [Ref.  25]  A  can  be  replaced  by  its 
definition,  yielding  an  ordinary  differential  equation  which 
admits  an  analytical  solution. 
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Spreiter  recently  extended  this  approach  to  the  oscillating 
airfoil  in  transonic  flow  [Ref.  28].   Start  with  the  oscillating 
flow  potential  equation  for  sonic  flow, 

K2 

x      Vx~s) 
iKx,o)  =  -  -1-   /  v(S) dS, 

■m  -l  -{x^b 

where , 

Q  =  A  +  2iK, 


X  =  constant,  replacing  the  variable  coefficient  (y+1)<J>vv 

AA 

With  ip(X,0),  calculate  ty   (X,0),  replace  X  by  the  function 
it  originally  represented,  and  solve  the  resulting  second 
order  differential  equation  for  i/;(X,0)  and  ^V(X,0). 

It  remains  to  be  investigated  whether  this  approach  could 
be  further  extended  to  study  nonsteady  interference  in  trans- 
onic flows.   It  appears  from  the  outset,  however,  that  the 
Collocation  method  is  suited  to  a  form  of  'local  linearization' 
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VII.   CONCLUSIONS 

The  unsteady,  sonic,  small  disturbance  potential  equation 
was  approximated,  following  the  approach  of  Teipel  [Ref.  30] 
and  Hosokawa  [Ref.  11],  by  replacing  the  nonconstant  coeffi- 
cients with  a  constant,  the  Oswatitsch  parameter.   The  re- 
sulting parabolic  equation  was  solved  by  Laplace  and  Fourier 
Transform  and  Collocation  techniques. 

Comparison  of  solutions  to  numerous  sonic  cascade  problems 
revealed  that  the  Transform  solutions  were  reliable  and  con- 
venient in  the  Fourier  Transform  case.   The  Collocation  tech- 
nique was  usable  with  qualifications,  growing  out  of  the  series 
approximation  and  numerical  integration  methods  utilized, 
which  degraded  its  accuracy  for  certain  parameters. 

Sonic  cascade  stability  boundaries  were  computed  utilizing 
the  Fourier  Transform  solution  and  interference  was  found  to 
be  strongly  destabilizing.   Such  results  must  be  qualified  by 
the  realization  that  they  were  obtained  by  the  solution  of  an 
approximated  equation,  though  the  trends  are  not  expected  to 
be  materially  affected. 

The  Collocation  method,  once  some  confidence  was  obtained 
from  its  comparison  with  the  other  solution  as  mentioned 
above,  was  applied  to  the  sonic  porous  wall  wind  tunnel  inter- 
ference problem.   Though  few  experimental  results  were  avail- 
able for  comparison  and  the  Drake  [Ref.  6]  approach  was 
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questioned,  the  trends  of  the  results  for  different  values 
of  porosity  were  as  expected  and  solid  wall  wind  tunnel  and 
sonic  free  jet  results  were  easily  obtainable  with  reasonable 
porosity  values. 

The  success  of  the  Collocation  solution  in  handling  a 
range  of  sonic  interference  problems  leads  to  the  considera- 
tion of  improvements  beyond  better  series  approximations  and 
numerical  techniques.   It  is  reasonable  to  expect  application 
of  the  Spreiter  [Ref.  28]  'local  linearization'  to  improve 
accuracy.   Further,  the  method  appears  to  be  capable  of  uti- 
lizing a  streamwise  variable  porosity  in  a  manner  like  that 
proposed  by  Sears  [Ref.  25]. 
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APPENDIX  A 
FORMULATION  OF  THE  LAPLACE  SOLUTION  ROUTINE 

A-I   Main  Program 

The  Main  routine  serves  as  the  controller  for  calculation 
of  potential  utilizing  the  Laplace  Transform  solution  derived 
in  Section  III-A.  % 

2 
K  (x-S)  _Q  (2np+p)2 

Tp(X,0  )  =  J —  E  ivn(S)e   [2e  ] 

-VttD?  -1     -VX^F   n=0 

Q  (2np+2p)2      Q  (2np)2 


"  V0 


(S)[e  ?   (X~S)   "  +  e  T   (X"S)]  }dS.      (111-2*0 


Calculations  are  made  for  controlling  variables  DK ,  DL, 
RHO,  and  SIGMA  (corresponding  to  K,  A,  p,  and  o  in  equation 
(111-24))  at  each  of  L  selected  points  along  the  chord. 

To  reduce  computational  redundancy,  the  Main  routine 
calculates  the  parameters  A,  C,  and  D  as  described  in  Section 
III-D.   Also,  the  boundary  condition  expressions 


v0(S)e10 


and 


v0(S) 


109 


are  treated  as  they  are  in  the  Collocation  routine  (Appendix 
C) 


vo(s)  ■  eoo +  eorx 


=    -1   -    iKX    ,  (A-I-l) 


eiav0(s)   =  e10  +  en.x  CA-I-2) 


The  Main  routine  reduces  equation  (111-22)  to 


POT(J)  =  R  *  I,  (A-I-3) 


where 


ri  =  —  ,     Q  =  X   +  2iK  ,  (A-I-4) 

■fnQ 


and 
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2         n    ,n_       ^2 


00 


i  *  i   {2.e10 


XExpC^  (X-S)  -  £^f+-4-] 


o"u  ;  "  ^      Cx^sj 


dS 


n=0    xu  .i         -tf^S 


X   Exp[^(X-S)  -  g  (2?fgl2]  dS 

+  2-eni   /  s  9 5 (x"s) 

11  -i         ^x^ 

X  Exp[g(X-S)  -  |  (2?f  g?)2] 

-  e  r  /  — -9 5 ObSj dS 

00  -l        -^x^ 

x  ^ynrx2rx  s^   Q  (2np)2i 

+    / « 2 u  b;   dS  ] 

-1        -/x^s1 

X   Exp[^(X-S)  -  f  i2njt|PL] 

-  en.  [  /  s  S 2 ^=^ —  ds 

01  -1  -/x^ 

x   exP[^(x-s)  -  §%4y] 

+   /  s  2 i! — yk^L_  dS]},  (a-i-5) 

-l        -Vx^s1 


which  is  calculated  by  subroutine  ISUM. 

The  Main  routine  uses  the  common  named  SUM  to  relay 
boundary  condition  and  integral  sum  results  between  itself 
and  subroutine  ISUM.   Common  HOLD  relays  parameter  information 
between  Main  and  subroutines  ISUM,  ZSUM,  and  INTCOR. 
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A-II   Subroutine  ISUM 

ISUM  is  called  by  Main  for  each  X  position  on  the  blade 
at  which  the  potential  is  to  be  evaluated.   It  sequentially 
calculates  each  term  in  the  summation  comprising  I  in  equation 
(A-I-5),  in  the  'n  loop'. 

ISUM  denotes  the  integrals  in  equation  (A-I-5)  as  the 
variables  Iji  where  j  =  1,2  and  i  =  0,1,2.  Since  ZSUM  is 
used  for  the  evaluation  of  the  integrals,  j  and  a  constant 

denote  the  power  of  U  required  in  the  integrand  (note  that 

i-1  5 
this  routine's  ZSUM  calculates  the  integrals  with  U 

1  -2  5 
vice  UJ      as  illustrated  in  Section  III-D) . 

Because  each  integral  differs  not  only  in  j  >   hut  in  the 

value  of  the  natural  exponent  values,  the  exponent  is 

considered  to  be  of  the  form 

Exp[K_(X-S)  -  £  ^1  , 

where  R0  =  2np,  Rl  =  2np+p,  R2  =  2np+2p.  Thus,  once  in  the 
'n  loop',  ISUM  calculates  R0 ,  Rl ,  and  R2  and  the  associated 
Bi  and  Ei  required  for  the  subsequent  integrations. 

The  integrations  are  conducted  sequentially  on  i.  Since 
Ij2  for  n-1  and  IjO  for  n  do  not  differ,  integration  for  IjO 
can  be  avoided  for  n  >  0. 

In  the  event  that  A  =  0,  correction  integrals  must  be 
calculated  as  outlined  in  Section  III-D.   Subroutine  INTCOR 
is  called  for  the  calculations  which  are  stored  matrix 
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COR(i,j,k).   The  corrections  are  to  be  calculated  but  once 
for  each  n;  thus,  NMAX  is  set  on  the  first  call  of  ISUM  and 
subsequent  calls  use  it  to  determine  whether  further  correc- 
tion calculations  are  required  (resetting  NMAX  higher). 

Accuracy  of  the  correction  integrals  is  dependent  upon 
the  magnitude  of  the  ratio  (E/U)/DU,  where  U  takes  on  the 
maximum  value  used  in  the  correction  integral,  i.e.,  U = XA 
(see  Section  III-D) .   If  we  desire  (E/U)/DU  >  10,000,  since 


E  =  |  (2np  +  kp)2,    k  =  0,1,2, 


then,  the  maximum  value  XA  can  be 

XA  =  (2n  +  k)p/100. 

Thus,  XA  is  calculated  for  each  n  and  i.   If  the  calculated 

XA  is  less  than  the  required  X0,  it  is  used.   If  it  is  greater, 

an  XA  is  chosen  smaller  than  X0. 

Of  course,  for  X   j-    0  and  R0  =  0.0,  XA  is  chosen  as 
1.0D-03  and  ZSUM  calculates  a  correction  (see  Section  III-D). 
Subsequent  integrations  utilize  the  arbitrarily  small 
XA  =  1.0D-08  as  a  compromise  between  computation  time  and 
accuracy. 

Convergence  of  the  series  for  I  is  signaled  when  successive 
terms  alter  it  by  less  than  DEL  x  100  percent  and  control  is 
returned  to  Main  with  the  value  for  I. 
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A-III.   Subroutine  INTCOR 

Called  by  ISUM  in  the  case  X  =  0  INTCOR  follows  the 
organization  of  Elder  [Ref.  7]  in  calling  the  corrections 
required  for  the  integrals  from  the  subroutine  COR. 

A-IV.   Subroutines  COR  and  ZSUM 

COR  is  a  simplification  of  the  more  general  correction 
routine  developed  for  the  Porous  Tunnel  Wall  Solution.   ZSUM 
is  an  adaptation  of  Elder's  [Ref.  7]  pyramided  Simpson's 
Rule  integration  routine  with  P  =  0  corrections  added  (see 
Section  III-D) . 

A-V.   Subroutine  CS 

This  is  a  library  routine  utilized  to  provide  Fresnel 
Integral  values  required  in  X   =  0  correction  integral 
calculations . 
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APPENDIX  B 
FORMULATION  OF  THE  FOURIER  SOLUTION  ROUTINE 

B-I   Main  Program 

The  Main  routine  satisfies  equation  (111-47)  derived  in 
Section  III-B  for  potential  at  selected  locations,  X,  along 
the  unit  blade  chord  and  equations  (111-49,  51,  and  55)  for 
pressure  coefficient,  total  lift  and  moment  coefficients. 

1|>(X,0)  =  -{(-1-iK  A  )[D.T1(X)+B-S1(X)]-iK1[D-T2(X)+E-S2(X)]>  , 

(111-47) 

Cp(X)  =  -2[^(X,0)  +  iK1tJ;(X,0)]  ,  (111-49) 

CL  =  2[>(1,0)  +  iK^Cl.O)]  ,  (III-5D 

CM  =  »2{iK1[F2(l,0)-F3(l,0)]+[^(l,0)-F2(l,0)]}-A1-CL, 

(IH-55) 

where  T„  and  S„  are  provided  by  Subroutine  SUM  as  described 
in  the  next  section  and  FN  are  successive  integrations  of  ty. 

Main  controls  the  calculation  of  TN  and  SN  with  Subroutine 
SUM  in  peculiar  manner.   Since   SN  and  TN  are  actually  infinite 
series,  Main  calls  SUM  to  produce  successive  terms  in  those 
series,  the  sums  of  which  are  S^  and  TN . 

Convergence  is  signaled  for  the  loop  of  index  n  in  Main 
when  successive  calculations  of  the  potential,  for  a  particular 
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chordwise  position,  differ  by  less  than  criterion  EPS. 

When  the  lift  and  moment  are  calculated,  the  n  loop  is 
continued,  if  necessary,  to  provide  the  successive  calcula- 
tions of  the  magnitude  and  phase  angle  of  the  moment  coeffi- 
cient until  convergence  is  detected  in  them. 

Due  to  the  difference  in  nondimensionalizations  between 
the  Fourier  and  other  solutions,  RKO,  DK ,  and  DL  must  be 
scaled,  as  shown  in  the  output,  to  conform  to  the  solution. 
Further,  the  results  for  the  potential  obtained  will  be  one 
half  those  of  the  other  solutions. 
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B-II.   Subroutine  SUM 

SUM  is  called  by  Main,  as  explained,  to  calculate  the 
series  for  SN  and  TN<   Recall  that  S  (X)  and  T-L(X)  are  from 
equations  (111-66  and  67).   Main  provides  the  coefficients 
to  the  series  as  Z. 

Since 

X 
bN+l  "  QJ         N 

and  likewise  for  T„ ,  then  the  dummy  variables  U„  in  SUM 
provide  the  appropriate  series  terms  for  S^r  or  T„  for  the 
value  of  J  (J  =  2  for  T^)  sent  to  it  in  the  calling  statement 

Logic  is  provided  to  properly  calculate  the  first  term 
of  the  series  for  TV,  for  the  case  K  =  0. 
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APPENDIX  C 
FORMULATION  OF  THE  COLLOCATION  SOLUTION  ROUTINE 

C-I.   Main  Program 

Because  the  derivations  have  certain  parallels  in  the 
integrals  to  be  evaluated,  much  of  this  routine  and  the 
Laplace  Transform  routine  are  similar. 

Main  calculates  the  necessary  variables  required  by  the 
integrating  routine,  the  boundary  condition  for  coefficient 
and  potential  generation,  and  the  step  size  between  colloca- 
tion points  (XINT).   With  that  done,  Subroutine  INTEGR  is 
called  to  calculate  all  integrals  required  in  subsequent 
steps.   Main  obtains  the  coefficients  from  Subroutine  COEFF, 
then,  implements  equation  (111-88)  and  calculates  the  potential 

Main's  organization  of  the  solution  for  the  potential  is 
also  applied  in  COEFF  and  is  due  to  the  series  chosen.   The 
J  loop  carries  the  calculations  down  the  chordwise  collocation 
points  and  the  I  loop  conducts  only  those  calculations  as 
allowed  by  the  chord  position  dependent  step  function  in  equa- 
tion (111-88).   For  any  X  position,  at  least  two  terms  must 
be  calculated,  the  n  =  0,  1  terms  of  the  series,  and  are  done 
when  1=1.   For  I  =  2,  no  calculations  need  be  made  and  it 
is  sidestepped  and  I  is  cycled  to  3. 
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C-II.   Subroutine  COEFF 

This  subroutine  is  called  by  Main  after  all  integrals 
have  been  calculated  for  the  purpose  of  solving  equations 
(111-86  and  87)  for  the  coefficients  of  the  series  representinj 
the  interference  source  distribution. 

For  n  =  0  and  1,  the  abovementioned  equations  can  be 
reduced  to  matrix  form, 


n=0 


n=l 


1 

X 

RIO 

RI1 

"v00 

"-BB0 

RIO 

RH 
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X 

v01 

-Bril 

1 

X 

RIO 

RI1 

V10 

-BBO 

_RI0 

RI1 

1 

x_ 

-vll. 

rBBl 

The  first  two  rows  are  calculated  at  the  first  collocation 
point,  selected  as  X  =  -1.0  +  XINT/2 . 0 ,  and  the  next  two  at 
the  second  point,  X  =  -1.0  +  1.5*XINT. 

The  above  matrix  equation  is  solved  by  a  complex  matrix 
inversion  routine  (CMINV)  and  a  multiplication  to  yield  the 

coefficients  NU(I,J),  I, J  =  1,2.  (i.e.,  vn(vvoi  ,vi  0 '  and  vll^* 
The  presence  of  the  step  function  in  the  equations  greatly 
simplifies  subsequent  calculations.   A  matrix  solution  is 
used  again  but  only  a  single  left  side  matrix  need  be  calcu- 
lated and  inverted.   As  successive  collocation  points  are 
reached,  all  previous  coefficients  have  been  calculated  and 
they  are  used  to  calculate  a  new  right  hand  side  vector. 
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C-III.   Subroutine  INTEGR 

Because  of  the  variety  of  integrals  required,  they  are 
calculated  as  soon  as  possible  and  stored  for  further  refer- 
ence in  the  matrix  INT(5,m)  carried  by  Common  INTG. 

Since  the  integrals  are  of  the  same  form  as  those  encoun- 
tered in  the  Laplace  solution,  the  same  solution  is  used,  with 
minor  differences  for  the  differing  exponent  J  range.   Sub- 
routine ZSUM,  borrowed  from  Elder  [Ref.  7  ]  and  modified  for 
B  =  E  =  0.0,  is  utilized  for  the  majority  of  the  integrations. 
In  the  event  that  the  Oswatitsch  parameter  is  taken  as  zero, 
Subroutine  COR  is  utilized  to  calculate  the  'correction 
integral ' . 

INTEGR  provides  five  different  integrals,  for  the  following 
sets  of  P  and  J  values, 
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for  the  ranges  0  to  X   =  -1.0  +  XINT/2.0  +  I      XINT  because 
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where 


U  =  X-S. 


Thence,  the  integrals  are  calculated  stepwise,  i.e.  the 
range  of  integration  starts  as  0  to  XQ,  then  Xn  to  X  ,  etc., 
until  all  necessary  integrals  are  obtained. 
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C-IV.   Subroutines  COR  and  ZSUM 

Subroutine  COR  is  an  adaptation  of  the  more  general 
routine  for  the  'correction  integral'  used  in  the  Porous 
Tunnel  Wall  program  described  in  Appendix  D-V. 

Subroutine  ZSUM  is  a  modification  of  Elder's  [Ref.  7] 
pyramided  Simpson's  Rule  integration  routine  to  include 
p  =  0  corrections  (see  Section  III-D) . 

C-V.   Subroutine  CMINV 

Subroutine  CMINV  inverts  the  complex  matrices  for 
Subroutine  COEFF.   It  is  an  outgrowth  of  a  Library  Subroutine 
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APPENDIX  D 
FORMULATION  OF  THE  SONIC  POROUS  WALL  WIND  TUNNEL  ROUTINE 

D-I.   Main  Program 

This  Main  program  is  very  similar  to  the  Main  program 
detailed  in  Appendix  C-I  for  the  Collocation  solution  of  the 
cascade  problem.   Provision  is  made  in  this  routine,  however, 
for  the  placement  of  the  'tunnel  wall'  at  Y  =  p/2  (which  will 
be  RHO  in  this  program).   Furthermore,  two  different  Oswatitsch 
parameters  may  be  utilized,  one  for  the  airfoil,  DL,  and  one 
for  the  wall,  DLW.   Because  the  wall  has  no  'isolated  blade' 
source  distribution,  TH10  and  TH11  do  not  appear.   Also,  de- 
rivatives are  made  available  so  as  to  calculate  the  pressure 
coefficient  as  well  as  the  velocity  potential.   Finally,  the 
variable,  SIGV/L  is  the  porosity  factor  (see  Section  V),  which 
can  be  made  variable  along  the  wall  (see  Subroutine  COEFF 
here) . 

The  Main  program  not  only  controls  the  calculation  of 
the  integrals  required,  by  calling  Subroutine  INTEGR,  and 
obtains  the  appropriate  interference  source  distribution 
series  coefficients  from  COEFF,  it  calculates  the  potential, 
POT,  and  the  pressure  coefficient,  CP,  at  the  selected 
locations  (collocation  points)  along  the  chord. 

From  Section  V  the  potential  on  the  airfoil  is 

iJ>(X,0)  =  ♦q(X.O)  +  ^(X.O)  +  (J)1(X,Y1=-p).         (D-l-1) 
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Also 


*X(X,0)    =    4>qX    (X,0)    +    *0X(X,0)    +    ^(X.Y^-p).  (D-I-2) 


The   pressure    coefficient    is 


Cp(X,0)    =    -2Df>    (X,0)    +    iK^(X,0)].  (D-I-3) 


By  introducing  the  step-wise-linear  series  used  in  the 
Collocation  solution  for  the  cascade  (equations  III-81  and 
82)  equation  V-2  becomes, 


*<x,t»  -  -V[0oo+voo]PIO°  +  [e01+v01]Pioi  +  J2V0nPioN} 

L-l 
+R3Cv10PI10+v11PIll   +      I    vlnPHN]    ,  (D-I-4) 


where, 

PI00    =         /  Exp[CldS,       PIoi    =    X      /X^£l£l  .        /XE2SP£C]dSj 

-1      -Vx^  -1    ix=s         -1     (x-s)  ±/d 


K2 

PION    =    PIOI   -    XN-PIOO,  C    =    £-(X-S), 

yo 

R2    =    p/2'^5      »  (D"I_5) 


15^ 


and,  for  PI10,  Pill,  and  PUN, 


K2        Ql   p2 
C  =  q-(X-S)  -  -jf  jj-gj 

R3  =  p/2  f±    • 


Then,  for  the  X  derivative  of  the  potential,  all  of  the 
PIxx  terms  of  the  potential  are  preceded  by  a  D,  denoting 
the  derivative.   The  appropriate  differentiated  integrals 

are  obtained  from  Section  V-A-3. 

2  2 

Y  Exp[5-(X-S)]      Exp[§-(X-S  )] 

DPIOO  =  ~   /  ^ dS  =  2 (D-I-6) 

dX  -1     -yx^  Y1^ 

etc  . 
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D-II.   Subroutine  COEFF 

As  the  Main  routine,  this  routine  differs  only  in  certain 
details  from  the  Collocation  routine  for  the  cascade.   Of 
greatest  importance  is  the  handling  of  the  wall  boundary 
condition  equation,  equation  (V-5) 


♦0Y(X»Y0=p)  +  WX»Y0=p)  +  *1Y(X'Y1=0)  + 


aw*{[*0X(X>Vp)  +  WX>Vp)  +  ♦lX^-Yl"0)] 


+  iK[^°(X,Y0=p)  +  (J)0(X,Y0-p)  +  4,1(X,Y1=0)]}  =  0.   (D-I-7) 


The  organization  of  the  solution  for  the  coefficients  is 
the  same  as  that  in  Appendix  C-II  save  for  the  inclusion  of 
the  SIGWL,  the  deletion  of  the  isolated  adjacent  blade  poten- 
tial and  normal  velocity  component,  and  a  slight  change  in 
notation . 

Illustrated  is  the  approach  for  varying  the  porosity. 
SIGWL  is  calculated  by  the  appropriate  function  in  COEFF  and 
used  as  necessary. 
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D-III.   Subroutine  INTEGR 

This  routine  is  similar  to  that  described  in  Appendix 
C-III  except  it  is  more  complex  to  calculate  the  additional 
integrals  required. 
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D-IV.   Subroutine  COR 

Subroutine  COR  calculates  the  'correction  integral' 
when  commanded  by  Subroutine  INTEGR,  using  equation  (111-99). 
The  'correction  integral'  is  obtained  by  summing  the  terms 
of  a  truncated  series,  here  SI(I)  and  CI(I),  I  =  0,1,2,  and  3 
SI(I)  and  CI(I)  are  integrals  of  the  form 


XA    /9 
/    Um/^  C?S(E/U)  dU 
0         Sln 


which  can  be  evaluated  by  integration  by  parts  and  use  of  the 
Fresnel  Integrals  (see  Section  III-D) .   Elder  [Ref.  7]  noted 
that  this  is  possible  and  reasonably  accurate  for  powers  of  U 
less  than  about  1.5.   In  any  event,  integrals  of  larger  powers 
are  very  small. 

COR  calculates  the  terms  by  liberal  use  of  equations 
(III-101  and  102).   It  must  be  noted  that  for  j  =  0  the  use  of 
equation  (III-101)  is  not  direct,  i.e.,  j  =  1  is  inserted  into 
(III-101)  and  the  resulting  integral  on  the  right  hand  side 
is  the  desired  j  =  0  integral 

XA 
J 
0 


XA 
;XA   y-5/2    sin(E/u)dU   =    ±1{    s     •  y-3/2    cos (£/u) dU+2XA~3/2cos (£/XA) 

cos  dru    ~.  sin  S-Ln 


Thus,  it  can  be  seen  why  the  involved  logic  is  contained  in 
COR  so  as  to  calculate  ocrrections  for  j  =  0  without  destroying 
the  basic  form  of  COR. 
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Of  course,  once  j  =  0  and  1  corrections  have  been 
calculated  no  further  calculations  of  terms  must  be  performed 
as  they  are  deemed  adequately  small  to  neglect. 
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D-V.   Subroutines  ZSUM,  CMINV,  and  CS 

These  subroutines  are  the  same  as  those  used  in  the 
Collocation  routine  for  the  cascade.   The  ZSUM  routine 
differs,  however,  in  that  the  terms  required  for  the 
'correction'  in  the  case  ?  =  0  are  precisely  equations 
(III-106),  to  provide  greater  flexibility. 
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